Bayesian GLM Part1

Author

Murray Logan

Published

06/07/2025

1 Preparations

Load the necessary libraries

library(tidyverse) # for data wrangling etc
library(rstanarm) # for fitting models in STAN
library(cmdstanr) # for cmdstan
library(brms) # for fitting models in STAN
library(standist) # for exploring distributions
library(coda) # for diagnostics
library(bayesplot) # for diagnostics
library(ggmcmc) # for MCMC diagnostics
library(DHARMa) # for residual diagnostics
library(rstan) # for interfacing with STAN
library(emmeans) # for marginal means etc
library(broom) # for tidying outputs
library(tidybayes) # for more tidying outputs
library(HDInterval) # for HPD intervals
library(ggeffects) # for partial plots
library(broom.mixed) # for summarising models
library(posterior) # for posterior draws
library(ggeffects) # for partial effects plots
library(patchwork) # for multi-panel figures
library(bayestestR) # for ROPE
library(see) # for some plots
library(easystats) # framework for stats, modelling and visualisation
library(INLA) # for approximate Bayes
library(INLAutils) # for additional INLA outputs
library(modelsummary) # for data and model summaries
theme_set(theme_grey()) # put the default ggplot theme back
source("helperFunctions.R")
bayesplot_theme_set(theme_bw(base_size = 8, base_family = "sans"))

2 Scenario

Here is an example from Fowler, Cohen, and Jarvis (1998). An agriculturalist was interested in the effects of fertilizer load on the yield of grass. Grass seed was sown uniformly over an area and different quantities of commercial fertilizer were applied to each of ten 1 m2 randomly located plots. Two months later the grass from each plot was harvested, dried and weighed. The data are in the file fertilizer.csv in the data folder.

Figure 1: Field of grass
Table 1: Format of the fertilizer.csv data file
FERTILIZER YIELD
25 84
50 80
75 90
100 154
125 148
... ...
Table 2: Description of the variables in the fertilizer data file
FERTILIZER: Mass of fertilizer (g.m-2) - Predictor variable
YIELD: Yield of grass (g.m-2) - Response variable

The aim of the analysis is to investigate the relationship between fertilizer concentration and grass yield.

3 Read in the data

We will start off by reading in the Fertilizer data. There are many functions in R that can read in a CSV file. We will use a the read_csv() function as it is part of the tidyverse ecosystem.

fert <- read_csv("../data/fertilizer.csv", trim_ws = TRUE)
Rows: 10 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (2): FERTILIZER, YIELD

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

After reading in a dataset, it is always a good idea to quickly explore a few summaries in order to ascertain whether the imported data are correctly transcribed. In particular, we should pay attention to whether there are any unexpected missing values and ensure that each variable (column) has the expected class (e.g. that variables we expected to be considered numbers are indeed listed as either <dbl> or <int> and not <char>).

glimpse(fert)
Rows: 10
Columns: 2
$ FERTILIZER <dbl> 25, 50, 75, 100, 125, 150, 175, 200, 225, 250
$ YIELD      <dbl> 84, 80, 90, 154, 148, 169, 206, 244, 212, 248
## Explore the first 6 rows of the data
head(fert)
# A tibble: 6 × 2
  FERTILIZER YIELD
       <dbl> <dbl>
1         25    84
2         50    80
3         75    90
4        100   154
5        125   148
6        150   169
str(fert)
spc_tbl_ [10 × 2] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ FERTILIZER: num [1:10] 25 50 75 100 125 150 175 200 225 250
 $ YIELD     : num [1:10] 84 80 90 154 148 169 206 244 212 248
 - attr(*, "spec")=
  .. cols(
  ..   FERTILIZER = col_double(),
  ..   YIELD = col_double()
  .. )
 - attr(*, "problems")=<externalptr> 
fert |> datawizard::data_codebook()
fert (10 rows and 2 variables, 2 shown)

ID | Name       | Type    | Missings |    Values |  N
---+------------+---------+----------+-----------+---
1  | FERTILIZER | numeric | 0 (0.0%) | [25, 250] | 10
---+------------+---------+----------+-----------+---
2  | YIELD      | numeric | 0 (0.0%) | [80, 248] | 10
-----------------------------------------------------
fert |> modelsummary::datasummary_skim()
Unique Missing Pct. Mean SD Min Median Max Histogram
FERTILIZER 10 0 137.5 75.7 25.0 137.5 250.0
YIELD 10 0 163.5 64.0 80.0 161.5 248.0

4 Exploratory data analysis

The individual responses (\(y_i\), observed yields) are each expected to have been independently drawn from normal (Gaussian) distributions (\(\mathcal{N}\)). These distributions represent all the possible yields we could have obtained at the specific (\(i^{th}\)) level of Fertilizer. Hence the \(i^{th}\) yield observation is expected to have been drawn from a normal distribution with a mean of \(\mu_i\).

Although each distribution is expected to come from populations that differ in their means, we assume that all of these distributions have the same variance (\(\sigma^2\)).

We need to supply priors for each of the parameters to be estimated (\(\beta_0\), \(\beta_1\) and \(\sigma\)). Whilst we want these priors to be sufficiently vague as to not influence the outcomes of the analysis (and thus be equivalent to the frequentist analysis), we do not want the priors to be so vague (wide) that they permit the MCMC sampler to drift off into parameter space that is both illogical as well as numerically awkward.

As a starting point, lets assign the following priors (how we arrive on these values will be outlined in sections on defining priors):

  • \(\beta_0\): Normal prior centred at 164 with a variance of 65
  • \(\beta_1\): Normal prior centred at 0 with a variance of 1
  • \(\sigma\): either a half cauchy centred at 0 with a variance of 2 or an exponential with parameter of 0.016

Note, when fitting models through either rstanarm or brms, the priors assume that the predictor(s) have been centred and are to be applied on the link scale. In this case the link scale is an identity.

Model formula:

\[ \begin{align} y_i &\sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i &= \beta_0 + \beta_1 x_i\\ \beta_0 &\sim{} \mathcal{N}(164,65)\\ \beta_1 &\sim{} \mathcal{N}(0,1)\\ \sigma &\sim{} \mathcal{t}(3,0,65)\\ OR\\ \sigma &\sim{} \mathcal{cauchy}(0,65)\\ OR\\ \sigma &\sim{} \mathcal{Exp}(0.016)\\ OR\\ \sigma &\sim{} \mathcal{gamma}(2,0.05)\\ \end{align} \]

We intend to explore the relationship between Yield and Fertiliser using simple linear regression. Such an analysis assumes:

  • Normality: each population is assumed to be normally distributed. This can be initially assessed by exploring the distribution of the response via either boxplots, violin plots or histograms (depending on sample size).
  • Homogeneity of variance: each population is assumed to be equally varied. This can be initially assessed by exploring the distribution of observed values around an imaginary line of best fit through the cloud of data formed by a scatterplot of Yield against Fertiliser. If the spread of points increases (or decreases) along the trend line, it is likely that the populations are not equally varied.
  • Linearity: in fitting a straight line through the cloud of data, we are assuming that a linear trend is appropriate. If the trend is not linear, then the inferred relationship may be miss-representative of the true relationship. In addition to an imaginary line, we could fit either a loess or linear smoother through the cloud to help us assess the likelihood of non-linearity.
  • Independence: to help ensure that the estimates are unbiased, we must assume that all observations are independent. In this case, we assume that the observations were collected from a randomised design in which the level of Fertiliser was applied randomly to the different patches of the field. If however, the researchers had simply moved along the field applying progressively higher Fertiliser concentrations, then there is a chance that they have introduced additional biases or confounding factors. Similarly, if the Fertiliser levels were applied in increasing doses over time, there might also be biases. In the absence of any spatial or temporal information associated with the collection of data, we cannot directly assess this assumption.
ggplot(fert, aes(y = YIELD, x = FERTILIZER)) +
  geom_point() +
  geom_smooth()

ggplot(fert, aes(y = YIELD, x = FERTILIZER)) +
  geom_point() +
  geom_smooth(method = "lm")

ggplot(fert, aes(y = YIELD)) +
  geom_boxplot(aes(x = 1))

ggplot(fert, aes(y = YIELD)) +
  geom_violin(aes(x = 1))

ggplot(fert, aes(x = YIELD)) +
  geom_histogram()

ggplot(fert, aes(x = YIELD)) +
  geom_density()

Conclusions:

  • there is no evidence of non-normality, non-homogeneity of variance or non-linearity
  • we have no initial reason to suspect that the assumptions will not be satisfied and thus it is reasonable to assume that the results will be reliable.

5 Fit the model

Note, for routines that take more than a couple of seconds to perform (such as most Bayesian models), it is a good idea to cache the Rmarkdown chunks in which the routine is performed. That way, the routine will only be run the first time and any objects generated will be stored for future use. Thereafter, provided the code has not changed, the routine will not be re-run. Rather, knitr will just retrieve the cached objects and continue on.

One of the most difficult aspects of performing Bayesian analyses is the specification of priors. For instances where there are some previous knowledge available and a desire to incorporate those data, the difficulty is in how to ensure that the information is incorporated correctly. However, for instances where there are no previous relevant information and so a desire to have the posteriors driven entirely by the new data, the difficulty is in how to define priors that are both vague enough (not bias results in their direction) and yet not so vague as to allow the MCMC sampler to drift off into unsupported regions (and thus get stuck and yield spurious estimates).

For early implementations of MCMC sampling routines (such as Metropolis Hasting and Gibbs), it was fairly common to see very vague priors being defined. For example, the priors on effects, were typically normal priors with mean of 0 and variance of 1e+06 (1,00,000). These are very vague priors. Yet for some samplers (e.g. NUTS), such vague priors can encourage poor behaviour of the sampler - particularly if the posterior is complex. It is now generally advised that priors should (where possible) be somewhat weakly informative and to some extent, represent the bounds of what are feasible and sensible estimates.

The degree to which priors influence an outcome (whether by having a pulling effect on the estimates or by encouraging the sampler to drift off into unsupported regions of the posterior) is dependent on:

  • the relative sparsity of the data - the larger the data, the less weight the priors have and thus less influence they exert.
  • the complexity of the model (and thus posterior) - the more parameters, the more sensitive the sampler is to the priors.

The sampled posterior is the product of both the likelihood and the prior - all of which are multidimensional. For most applications, it would be vertically impossible to define a sensible multidimensional prior. Hence, our only option is to define priors on individual parameters (e.g. the intercept, slope(s), variance etc) and to hope that if they are individually sensible, they will remain collectively sensible.

So having (hopefully) impressed upon the notion that priors are an important consideration, I will now attempt to synthesise some of the approaches that can be employed to arrive at weakly informative priors that have been gleaned from various sources. Largely, this advice has come from the following resources:

  • https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations
  • http://svmiller.com/blog/2021/02/thinking-about-your-priors-bayesian-analysis/

I will outline some of the current main recommendations before summarising some approaches in a table.

  • weakly informative priors should contain enough information so as to regularise (discourage unreasonable parameter estimates whilst allowing all reasonable estimates).
  • for effects parameters on scaled data, an argument could be made for a normal distribution with a standard deviation of 1 (e.g. normal(0,1)), although klsome prefer a t distribution with 3 degrees of freedom and standard deviation of 1 (e.g. student_t(3,0,1)) - apparently a flatter t is a more robust prior than a normal as an uninformative prior…
  • for unscaled data, the above priors can be scaled by using the standard deviation of the data as the prior standard deviation (e.g. student_t(3,0,sd(y)), or sudent_t(3,0,sd(y)/sd(x)))
  • for priors of hierachical standard deviations, priors should encourage shrinkage towards 0 (particularly if the number of groups is small, since otherwise, the sampler will tend to be more responsive to “noise”).
Family Parameter brms rstanarm
Gaussian Intercept student_t(3,median(y),mad(y)) normal(mean(y),2.5*sd(y))
‘Population effects’ (slopes, betas) flat, improper priors normal(0,2.5*sd(y)/sd(x))
Sigma student_t(3,0,mad(y)) exponential(1/sd(y))
‘Group-level effects’ student_t(3,0,mad(y)) decov(1,1,1,1)
Correlation on group-level effects ljk_corr_cholesky(1)
Poisson Intercept student_t(3,median(y),mad(y)) normal(mean(y),2.5*sd(y))
‘Population effects’ (slopes, betas) flat, improper priors normal(0,2.5*sd(y)/sd(x))
‘Group-level effects’ student_t(3,0,mad(y)) decov(1,1,1,1)
Correlation on group-level effects ljk_corr_cholesky(1)
Negative binomial Intercept student_t(3,median(y),mad(y)) normal(mean(y),2.5*sd(y))
‘Population effects’ (slopes, betas) flat, improper priors normal(0,2.5*sd(y)/sd(x))
Shape gamma(0.01, 0.01) exponential(1/sd(y))
‘Group-level effects’ student_t(3,0,mad(y)) decov(1,1,1,1)
Correlation on group-level effects ljk_corr_cholesky(1)

Notes:

brms

https://github.com/paul-buerkner/brms/blob/c2b24475d727c8afd8bfc95947c18793b8ce2892/R/priors.R

  1. In the above, for non-Gaussian families, y is first transformed according to the family link. If the family link is log, then 0.1 is first added to 0 values.
  2. in brms the minimum standard deviation for the Intercept prior is 2.5
  3. in brms the minimum standard deviation for group-level priors is 10.

rstanarm

http://mc-stan.org/rstanarm/articles/priors.html

  1. in rstanarm priors on standard deviation and correlation associated with group-level effects are packaged up into a single prior (decov which is a decomposition of the variance and covariance matrix).

For the purpose of comparing the frequentist and Bayesian outcomes, it might be useful to start by fitting the frequentist simple linear model. We will also fit this model with the predictor (fertilizer concentration) centred.

summary(fert.lm <- lm(YIELD ~ FERTILIZER, data = fert))

Call:
lm(formula = YIELD ~ FERTILIZER, data = fert)

Residuals:
   Min     1Q Median     3Q    Max 
-22.79 -11.07  -5.00  12.00  29.79 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 51.93333   12.97904   4.001  0.00394 ** 
FERTILIZER   0.81139    0.08367   9.697 1.07e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 19 on 8 degrees of freedom
Multiple R-squared:  0.9216,    Adjusted R-squared:  0.9118 
F-statistic: 94.04 on 1 and 8 DF,  p-value: 1.067e-05
sigma(fert.lm)
[1] 18.99936
summary(fert.lm <- lm(YIELD ~ center(FERTILIZER), data = fert))

Call:
lm(formula = YIELD ~ center(FERTILIZER), data = fert)

Residuals:
   Min     1Q Median     3Q    Max 
-22.79 -11.07  -5.00  12.00  29.79 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        163.50000    6.00813  27.213 3.58e-09 ***
center(FERTILIZER)   0.81139    0.08367   9.697 1.07e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 19 on 8 degrees of freedom
Multiple R-squared:  0.9216,    Adjusted R-squared:  0.9118 
F-statistic: 94.04 on 1 and 8 DF,  p-value: 1.067e-05
sigma(fert.lm)
[1] 18.99936

In rstanarm, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over-fitting) and help stabilise the computations.

fert.rstanarm <- stan_glm(YIELD ~ FERTILIZER,
  data = fert,
  iter = 5000,
  warmup = 1000,
  chains = 3,
  thin = 5,
  refresh = 0
)

In the above:

  • the formula and data arguments should be familiar as they are the same as for many models in R
  • iter: indicates the number of MCMC iterations to perform per chain
  • warmup: indicates how many of the initial MCMC samples to exclude. During the warmup stage, the MCMC sampler is tuning its sampling operations to allow it to determine the best way to create the chain. While it is in this discovery phase, the estimates it produces might be unrepresentative.
  • chains: indicates the number of independent chains to run. The more chains, the more random starts. The point of multiple random starts is to minimise the potential for the starting point to dictate the outcome. It is possible that a chain might get stuck sampling a single feature of the posterior and not discover additional features. Additional random starting points should hopefully alleviate this. Either 3 or 4 is typical.
  • thin: indicates the rate at which the MCMC samples should be thinned in order to eliminate auto-correlation between MCMC samples. When the step lengths between MCMC samples are small, the resulting parameter estimates will be auto-correlated (since samples close by in the chain will be more similar than those far apart). The rate of thinning necessary will depend on the degree of auto-correlation. For stan models, a thinning rate of 5 is a good starting point.
  • refresh: indicates how often the terminal display should be updated. When running a large model, this can be useful as it provides a form of progress. However, when the routine is within an Rmarkdown document, it just results in the output including each line of progress and a lot of space is taken up unnecessarily

Having allowed rstanarm to formulate its own weakly informative priors, it is a good idea to explore what they are. Firstly, out of curiosity, it might be interesting to see what it has chosen. However, more importantly, we need to be able to document what the priors were and the rstanarm development team make it very clear that there is no guarantee that the default priors will remain the same into the future.

prior_summary(fert.rstanarm)
Priors for model 'fert.rstanarm' 
------
Intercept (after predictors centered)
  Specified prior:
    ~ normal(location = 164, scale = 2.5)
  Adjusted prior:
    ~ normal(location = 164, scale = 160)

Coefficients
  Specified prior:
    ~ normal(location = 0, scale = 2.5)
  Adjusted prior:
    ~ normal(location = 0, scale = 2.1)

Auxiliary (sigma)
  Specified prior:
    ~ exponential(rate = 1)
  Adjusted prior:
    ~ exponential(rate = 0.016)
------
See help('prior_summary.stanreg') for more details

Note, consistent with good practice, rstanarm will centre all categorical predictors (whether the user centres them or not). Whilst this does not impact on estimated slope(s), it does alter the intercept (arguably for the better). Since all priors are applied on the scale of the parameters in the linear predictor (e.g. the intercept and slopes), the prior on intercept must be defined as if the predictors have been centred.

Nevertheless, rstanarm will return a posterior for intercept that is consistent with how the user had defined the linear predictor - that is, if the user does not centre the predictors, then rstanarm will return the intercept posterior on a scale consistent with as if the predictors had not been centred, even though they were automatically scaled internally.

When defining the spread of priors, rstanarm prefers to do so relative to the scaling of the data. For example, for the slope, it takes the ratio of the standard deviations of the response and the predictor and multiplies them by a factor of 2.5 (by default). Hence it is possible to define priors either in relative terms (by providing a scaling factor) or in absolute terms, by providing the desired priors along with the autoscale = FALSE argument.

The above prior summary tells us:

  • for the intercept, it is using a normal prior with a mean of 164 and a standard deviation of 2.5. This mean value is based on the mean of the response variable (rounded up). The 2.5 is used for all intercepts. The scale of the prior is then adjusted proportional to the variability in the response, by multiplying 2.5 by the standard deviation of the response (and rounding up).
mean(fert$YIELD)
[1] 163.5
sd(fert$YIELD)
[1] 63.97439
2.5 * sd(fert$YIELD)
[1] 159.936
  • for the coefficients (in this case, just the slope), the default prior is a normal prior centred around 0 with a standard deviation of 2.5. This is then adjusted for the scale of the data by multiplying the 2.5 by the ratio of the standard deviation of the response to the standard deviation of the predictor (then rounded).
2.5 * sd(fert$YIELD) / sd(fert$FERTILIZER)
[1] 2.113004
  • the auxiliary prior represents whatever additional prior is required for the nominated model. In the case of a Gaussian model, this will be \(\sigma\), for negative binomial, it will be the reciprocal of dispersion, for gamma, it will be shape, etc . By default in rstanarm, this is a exponential with a rate of 1 which is then adjusted by division with the standard deviation of the response.
1 / sd(fert$YIELD)
[1] 0.01563126

One way to assess the priors is to have the MCMC sampler sample purely from the prior predictive distribution without conditioning on the observed data. Doing so provides a glimpse at the range of predictions possible under the priors. On the one hand, wide ranging predictions would ensure that the priors are unlikely to influence the actual predictions once they are conditioned on the data. On the other hand, if they are too wide, the sampler is being permitted to traverse into regions of parameter space that are not logically possible in the context of the actual underlying ecological context. Not only could this mean that illogical parameter estimates are possible, when the sampler is traversing regions of parameter space that are not supported by the actual data, the sampler can become unstable and have difficulty.

We can draw from the prior predictive distribution instead of conditioning on the response, by updating the model and indicating prior_PD=TRUE. After refitting the model in this way, we can plot the predictions to gain insights into the range of predictions supported by the priors alone.

fert.rstanarm1 <- update(fert.rstanarm, prior_PD = TRUE)
ggpredict(fert.rstanarm1) |> plot(show_data = TRUE)
Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

Conclusions:

  • we see that the range of predictions is fairly wide and the slope could range from strongly negative to strongly positive.
  • at the same time, massive yields are not supported, so there is some constraints on the parameters space
  • nevertheless, negative yields are considered possible, even though these are clearly not logical.

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):

  • \(\beta_0\): normal centred at 164 with a standard deviation of 65
  • \(\beta_1\): normal centred at 0 with a standard deviation of 2
  • \(\sigma\): half-t (df = 3) centred at 0 with a standard deviation of 65

We can visualise each of the above as:

standist::visualize("normal(164, 65)", xlim = c(0, 300)) +
  standist::visualize("normal(0, 2)", xlim = c(-10, 10))

standist::visualize("student_t(3, 0, 65)",
  "gamma(2,1)",
  "exponential(0.016)",
  xlim = c(-10, 50)
)

standist::visualize("student_t(3, 0, 65)",
  "exponential(0.016)",
  xlim = c(-10, 300)
)

I will also overlay the raw data for comparison.

fert.rstanarm2 <- stan_glm(YIELD ~ FERTILIZER,
  data = fert,
  prior_intercept = normal(164, 65, autoscale = FALSE),
  prior = normal(0, 2, autoscale = FALSE),
  prior_aux = student_t(3, 0, 65, autoscale = FALSE),
  prior_PD = TRUE,
  iter = 5000, warmup = 1000,
  chains = 3, cores = 3,
  thin = 5, refresh = 0
)
ggpredict(fert.rstanarm2) |>
  plot(show_data = TRUE)
Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

Now lets refit, conditioning on the data.

fert.rstanarm3 <- update(fert.rstanarm2, prior_PD = FALSE)
posterior_vs_prior(fert.rstanarm3,
  color_by = "vs", group_by = TRUE,
  facet_args = list(scales = "free_y")
)

Drawing from prior...

Conclusions:

  • in each case, the prior is substantially wider than the posterior, suggesting that the posterior is not biased towards the prior.
ggpredict(fert.rstanarm3) |> plot(show_data = TRUE)
Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over-fitting) and help stabilise the computations.

Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.

fert.brm <- brm(bf(YIELD ~ FERTILIZER),
  data = fert,
  iter = 5000,
  warmup = 1000,
  chains = 3,
  thin = 5,
  refresh = 0,
  backend = "rstan"
)
Compiling Stan program...
Start sampling

In the above:

  • the formula and data arguments should be familiar as they are the same as for many models in R. That said, the formula specification of brms can define models that are not possible by most other routines. To facilitate this enhanced functionality, we usually define a brms formula within its own bf() function along with the family (in this case, it is Gaussian, which is the default and therefore can be omitted.)
  • iter: indicates the number of MCMC iterations to perform per chain
  • warmup: indicates how many of the initial MCMC samples to exclude. During the warmup stage, the MCMC sampler is tuning its sampling operations to allow it to determine the best way to create the chain. While it is in this discovery phase, the estimates it produces might be unrepresentative.
  • chains: indicates the number of independent chains to run. The more chains, the more random starts. The point of multiple random starts is to minimise the potential for the starting point to dictate the outcome. It is possible that a chain might get stuck sampling a single feature of the posterior and not discover additional features. Additional random starting points should hopefully alleviate this. Either 3 or 4 is typical.
  • thin: indicates the rate at which the MCMC samples should be thinned in order to eliminate auto-correlation between MCMC samples. When the step lengths between MCMC samples are small, the resulting parameter estimates will be auto-correlated (since samples close by in the chain will be more similar than those far apart). The rate of thinning necessary will depend on the degree of auto-correlation. For stan models, a thinning rate of 5 is a good starting point.
  • refresh: indicates how often the terminal display should be updated. When running a large model, this can be useful as it provides a form of progress. However, when the routine is within an Rmarkdown document, it just results in the output including each line of progress and a lot of space is taken up unnecessarily

The returned object (a list) contains a range of objects. The output from rstan is contained in the fit list item.

fert.brm |> names()
 [1] "formula"   "data"      "prior"     "data2"     "stanvars"  "model"    
 [7] "save_pars" "algorithm" "backend"   "threads"   "opencl"    "stan_args"
[13] "fit"       "basis"     "criteria"  "file"      "version"   "family"   
[19] "autocor"   "ranef"     "cov_ranef" "stan_funs" "data.name"

Having allowed brms to formulate its own weakly informative priors, it is a good idea to explore what they are. Firstly, out of curiosity, it might be interesting to see what it has chosen. However, more importantly, we need to be able to document what the priors were and the brms development team make it very clear that there is no guarantee that the default priors will remain the same into the future.

prior_summary(fert.brm)
                     prior     class       coef group resp dpar nlpar lb ub       source
                    (flat)         b                                             default
                    (flat)         b FERTILIZER                             (vectorized)
 student_t(3, 161.5, 90.4) Intercept                                             default
     student_t(3, 0, 90.4)     sigma                                   0         default

This tells us:

  • for the intercept, it is using a student t (flatter normal) prior with a mean of 161.5 and a standard deviation of 90.4. These mean and standard deviation values are based on the median and median absolute deviation of the response variable (rounded).

Note, all priors are on the link scale and on centered predictiors

median(fert$YIELD)
[1] 161.5
mad(fert$YIELD)
[1] 90.4386

mad is the median absolute deviation. mad is to var as median is to mean.

standist::visualize("student_t(3,161.5,90.4)", xlim = c(-100, 1000))

  • for the beta coefficients (in this case, just the slope), the default prior is a improper flat prior. A flat prior essentially means that any value between negative infinity and positive infinity are equally likely. Whilst this might seem reckless, in practice, it seems to work reasonably well for non-intercept beta parameters.

  • the prior on sigma is also a student t (although only the positive half is used in the stan code), with mean of 0 and standard deviation of 90.4.

standist::visualize("student_t(3,0,90.4)", xlim = c(-10, 500))

One way to assess the priors is to have the MCMC sampler sample purely from the prior predictive distribution without conditioning on the observed data. Doing so provides a glimpse at the range of predictions possible under the priors. On the one hand, wide ranging predictions would ensure that the priors are unlikely to influence the actual predictions once they are conditioned on the data. On the other hand, if they are too wide, the sampler is being permitted to traverse into regions of parameter space that are not logically possible in the context of the actual underlying ecological context. Not only could this mean that illogical parameter estimates are possible, when the sampler is traversing regions of parameter space that are not supported by the actual data, the sampler can become unstable and have difficulty.

In brms, we can inform the sampler to draw from the prior predictive distribution instead of conditioning on the response, by running the model with the sample_prior = 'only' argument. Unfortunately, this cannot be applied when there are flat priors (since the posteriors will necessarily extend to negative and positive infinity). Therefore, in order to use this useful routine, we need to make sure that we have defined a proper prior for all parameters.

Lets try a Gaussian (normal) distribution for the slope. We will start with what is likely to be a very wide distribution (wide with respect to what we know about the likely slope). rstanarm sets the default for such population-level effects at 2.5 times the ratio of standard deviations of the response and predictor. The reason that 2.5 is used is that in a standard normal distribution, a standard deviation of 2.5 encapsulates approximately 99% of the values. Hence, a standard deviation of 2.5 would represent reasonably wide expectations (priors). Given that the input data are unlikely to be standardised, we can instead standardise the priors by multiplying the 2.5 by the ratio of standard deviations of the response and predictor.

We could adopt similar logic here, but rather than use standard deviation, we could use median absolute deviation to be consistent with other prior specifications in brms. In this case, the prior standard deviation would be 2.5 * mad(fert$YIELD)/mad(fert$FERTILIZER) = 2.44 (which I will round to 2.5).

standist::visualize("normal(0, 2.5)", "student_t(3, 0, 2.5)", xlim = c(-50, 50))

fert.brm1 <- brm(bf(YIELD ~ FERTILIZER),
  data = fert,
  prior = prior(student_t(3, 0, 2.5), class = "b"),
  sample_prior = "only",
  iter = 5000,
  warmup = 1000,
  chains = 3,
  thin = 5,
  backend = "rstan",
  refresh = 0
)
Compiling Stan program...
Start sampling
fert.brm1 |>
  ggpredict() |>
  plot(show_data = TRUE)
Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

fert.brm1 |>
  ggemmeans(~FERTILIZER) |>
  plot(show_data = TRUE)
Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

fert.brm1 |>
  conditional_effects() |>
  plot(points = TRUE)

Conclusions:

  • we see that the range of predictions is fairly wide and the slope could range from strongly negative to strongly positive.
  • at the same time, massive yields are not supported, so there is some constraints on the parameters space
  • nevertheless, negative yields are considered possible, even though these are clearly not logical.

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (which I have mainly plucked out of thin air):

  • \(\beta_0\): normal centred at 164 with a standard deviation of 65
    • mean of 164: since median(fert$YIELD)
    • sd pf 65: since mad(fert$YIELD)
  • \(\beta_1\): normal centred at 0 with a standard deviation of 2.5
    • sd of 2: since 2.5*(mad(fert$YIELD)/mad(fert$FERTILIZER))
  • \(\sigma\): half-t centred at 0 with a standard deviation of 65 OR
    • sd pf 65: since mad(fert$YIELD)
  • \(\sigma\): gamma with shape parameters of 2 and 1
standist::visualize("normal(164, 65)", xlim = c(-100, 300))

standist::visualize("normal(0, 90)", xlim = c(-10, 10))

standist::visualize("student_t(3, 0, 65)",
  "normal(0, 65)",
  "gamma(2,0.05)",
  xlim = c(0, 100)
)

I will also overlay the raw data for comparison.

form <- bf(YIELD ~ scale(FERTILIZER, scale = FALSE))
priors <- prior(normal(162, 90.4), class = "Intercept") +
  prior(normal(0, 2), class = "b") +
  prior(student_t(3, 0, 90.4), class = "sigma")
fert.brm2 <- brm(form,
  data = fert,
  prior = priors,
  sample_prior = "only",
  iter = 5000,
  warmup = 1000,
  chains = 3, cores = 3,
  thin = 5,
  backend = "rstan",
  refresh = 0
)
Compiling Stan program...
Start sampling
fert.brm2 |>
  ggpredict() |>
  plot(show_data = TRUE)
Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

fert.brm2 |>
  ggemmeans(~FERTILIZER) |>
  plot(show_data = TRUE)
Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

fert.brm2 |>
  conditional_effects() |>
  plot(points = TRUE)

fert.brm3 <- update(fert.brm2, sample_prior = "yes", refresh = 0)
The desired updates require recompiling the model
Compiling Stan program...
Start sampling
# OR
fert.brm3 <- brm(bf(YIELD ~ scale(FERTILIZER, scale = FALSE)),
  data = fert,
  prior = priors,
  sample_prior = "yes",
  iter = 5000,
  warmup = 1000,
  chains = 3, cores = 3,
  thin = 5,
  backend = "rstan",
  refresh = 0
)
Compiling Stan program...
Start sampling
fert.brm3 |>
  ggpredict() |>
  plot(show_data = TRUE)
Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

fert.brm3 |>
  ggemmeans(~FERTILIZER) |>
  plot(show_data = TRUE)
Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

fert.brm3 |>
  conditional_effects() |>
  plot(points = TRUE)

When we have indicated that the posterior should be informed by both the prior and the posterior, both prior (governed by priors alone) and posterior (governed by both priors and data/likelihood) draws are returned. These can be compared by exploring the probabilities associated with specific hypotheses - the most obvious of which is that of no effect (that the parameter = 0).

When doing so, ideally the posterior should be distinct from the prior. If this is not the case, it might suggest that the posteriors are being too strongly driven by the priors.

fert.brm3 |> get_variables()
 [1] "b_Intercept"                   "b_scaleFERTILIZERscaleEQFALSE"
 [3] "sigma"                         "Intercept"                    
 [5] "prior_Intercept"               "prior_b"                      
 [7] "prior_sigma"                   "lprior"                       
 [9] "lp__"                          "accept_stat__"                
[11] "stepsize__"                    "treedepth__"                  
[13] "n_leapfrog__"                  "divergent__"                  
[15] "energy__"                     
## fert.brm3 |> hypothesis('Intercept=0', class='b') |> plot
fert.brm3 |>
  hypothesis("scaleFERTILIZERscaleEQFALSE = 0") |>
  plot()

## fert.brm3 |> hypothesis('scaleFERTILIZERscaleEQFALSE=0') |> plot
fert.brm3 |>
  hypothesis("sigma = 0", class = "") |>
  plot()

Unfortunately, it is not possible to do this comparison sensibly for the intercept. The reason for this is that the prior for intercept was applied to an intercept that is associated with centred continuous predictors (predictors are centred behind the scenes). Since we did not centre the predictor, the intercept returned is as if uncentred. Hence, the prior and posterior are on different scales.

fert.brm3 |> get_variables()
 [1] "b_Intercept"                   "b_scaleFERTILIZERscaleEQFALSE"
 [3] "sigma"                         "Intercept"                    
 [5] "prior_Intercept"               "prior_b"                      
 [7] "prior_sigma"                   "lprior"                       
 [9] "lp__"                          "accept_stat__"                
[11] "stepsize__"                    "treedepth__"                  
[13] "n_leapfrog__"                  "divergent__"                  
[15] "energy__"                     
fert.brm3 |> SUYR_prior_and_posterior()

fert.brm3 |> standata()
$N
[1] 10

$Y
 [1]  84  80  90 154 148 169 206 244 212 248

$K
[1] 2

$Kc
[1] 1

$X
   Intercept scaleFERTILIZERscaleEQFALSE
1          1                      -112.5
2          1                       -87.5
3          1                       -62.5
4          1                       -37.5
5          1                       -12.5
6          1                        12.5
7          1                        37.5
8          1                        62.5
9          1                        87.5
10         1                       112.5
attr(,"assign")
[1] 0 1

$prior_only
[1] 0

attr(,"class")
[1] "standata" "list"    
fert.brm3 |> stancode()
// generated with brms 2.22.0
functions {
}
data {
  int<lower=1> N;  // total number of observations
  vector[N] Y;  // response variable
  int<lower=1> K;  // number of population-level effects
  matrix[N, K] X;  // population-level design matrix
  int<lower=1> Kc;  // number of population-level effects after centering
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
  matrix[N, Kc] Xc;  // centered version of X without an intercept
  vector[Kc] means_X;  // column means of X before centering
  for (i in 2:K) {
    means_X[i - 1] = mean(X[, i]);
    Xc[, i - 1] = X[, i] - means_X[i - 1];
  }
}
parameters {
  vector[Kc] b;  // regression coefficients
  real Intercept;  // temporary intercept for centered predictors
  real<lower=0> sigma;  // dispersion parameter
}
transformed parameters {
  real lprior = 0;  // prior contributions to the log posterior
  lprior += normal_lpdf(b | 0, 2);
  lprior += normal_lpdf(Intercept | 162, 90.4);
  lprior += student_t_lpdf(sigma | 3, 0, 90.4)
    - 1 * student_t_lccdf(0 | 3, 0, 90.4);
}
model {
  // likelihood including constants
  if (!prior_only) {
    target += normal_id_glm_lpdf(Y | Xc, Intercept, b, sigma);
  }
  // priors including constants
  target += lprior;
}
generated quantities {
  // actual population-level intercept
  real b_Intercept = Intercept - dot_product(means_X, b);
  // additionally sample draws from priors
  real prior_b = normal_rng(0,2);
  real prior_Intercept = normal_rng(162,90.4);
  real prior_sigma = student_t_rng(3,0,90.4);
  // use rejection sampling for truncated priors
  while (prior_sigma < 0) {
    prior_sigma = student_t_rng(3,0,90.4);
  }
}

In INLA, the default priors are designed to be diffuse or weak. They are chosen to provide moderate regularisation (to help prevent over-fitting) and help stabilise the computations.

** Note as of mid 2022, the INLA developers have intentionally disabled the calculation of marginals in order to make the models faster**. Unfortunately, this prevents the calculation of residuals from the INLAutils package. So I will turn them on… More info can be viewed here https://groups.google.com/g/r-inla-discussion-group/c/2E9z3QfmeXw/m/XB0PoLqEBAAJ?pli=1

fert.inla <- inla(YIELD ~ FERTILIZER,
  data = fert,
  family = "gaussian",
  control.predictor = list(compute = TRUE, link = 1),
  control.compute = list(
    config = TRUE,
    dic = TRUE, waic = TRUE, cpo = TRUE,
    return.marginals.predictor = TRUE
  )
)

In the above:

  • the formula, data and family arguments should be familiar as they are the same as for other models in R.

  • control.compute: allows us to indicate what additional actions should be performed during the model fitting. In this case, we have indicated:

    • dic: Deviance information criterion
    • waic: Wantanabe information creterion
    • cpo: out-of-sample estimates (measures of fit)
    • config: return the full configuration - to allow drawing from the posterior.
fert.inla |> names()
 [1] "names.fixed"                 "summary.fixed"              
 [3] "marginals.fixed"             "summary.lincomb"            
 [5] "marginals.lincomb"           "size.lincomb"               
 [7] "summary.lincomb.derived"     "marginals.lincomb.derived"  
 [9] "size.lincomb.derived"        "mlik"                       
[11] "cpo"                         "gcpo"                       
[13] "po"                          "waic"                       
[15] "residuals"                   "model.random"               
[17] "summary.random"              "marginals.random"           
[19] "size.random"                 "summary.linear.predictor"   
[21] "marginals.linear.predictor"  "summary.fitted.values"      
[23] "marginals.fitted.values"     "size.linear.predictor"      
[25] "summary.hyperpar"            "marginals.hyperpar"         
[27] "internal.summary.hyperpar"   "internal.marginals.hyperpar"
[29] "offset.linear.predictor"     "model.spde2.blc"            
[31] "summary.spde2.blc"           "marginals.spde2.blc"        
[33] "size.spde2.blc"              "model.spde3.blc"            
[35] "summary.spde3.blc"           "marginals.spde3.blc"        
[37] "size.spde3.blc"              "logfile"                    
[39] "misc"                        "dic"                        
[41] "mode"                        "joint.hyper"                
[43] "nhyper"                      "version"                    
[45] "Q"                           "graph"                      
[47] "ok"                          "cpu.intern"                 
[49] "cpu.used"                    "all.hyper"                  
[51] ".args"                       "call"                       
[53] "model.matrix"               

Having allowed INLA to formulate its own “minimally informative” priors, it is a good idea to explore what they are. Firstly, out of curiosity, it might be interesting to see what it has chosen. However, more importantly, we need to be able to document what the priors were and the INLA development team make it very clear that there is no guarantee that the default priors will remain the same into the future.

In calcutating the posterior mode of hyperparameters, it is efficient to maximising the sum of the (log)-likelihood and the (log)-prior, hence, priors are defined on a log-scale. The canonical prior for variance is the gamma prior, hence in INLA, this is a loggamma.

They are also defined according to their mean and precision (inverse-scale, rather than variance). Precision is \(1/\sigma\).

To explore the default priors used by INLA, we can issue the following on the fitted model:

inla.priors.used(fert.inla)
section=[family]
    tag=[INLA.Data1] component=[gaussian]
        theta1:
            parameter=[log precision]
            prior=[loggamma]
            param=[1e+00, 5e-05]
section=[fixed]
    tag=[(Intercept)] component=[(Intercept)]
        beta:
            parameter=[(Intercept)]
            prior=[normal]
            param=[0, 0]
    tag=[FERTILIZER] component=[FERTILIZER]
        beta:
            parameter=[FERTILIZER]
            prior=[normal]
            param=[0.000, 0.001]

The above indicates:

  • the prior on the overal family (log) precision (\(1/sigma\)) is a log-gamma with shape of 1 and scale of 0.00005. Since this is a log-gamma on log-precision, it is equivalent to a gamma (shape = \(1\), scale = \(10^-5\)) on precision. Such a distribution is very wide and has a mean of (\(shape/scale =\) 2^{4}) and variance of (\(shape/scale^2 =\) 4^{8}).
standist::visualize("gamma(1, 0.00005)", xlim=c(-100,100000))

  • the prior on the intercept is a Gaussian (normal) distribution with mean of 0 and precision of 0 - this is effectively a uniform distribution from \(-\infty\) to \(\infty\).
  • the prior on the slope is a Gaussian (normal) distribution with a mean of 0 and a precision of 0.001 (which equates to a standard deviation of \((1/0.001)^{0.5}\)=31.6227766. This is also a very wide distribution (consider what a slope large in magnitude of 100 means in the context of the rate of grass yield increase per unit increase in fertiliser concentration).
standist::visualize("normal(0, 31)", xlim=c(-100,100))

Family variance

In order to drive the specification of weekly informative priors for variance components (a prior that allows the data to dictate the parameter estimates whilst constraining the range of the marginal distribution within sensible bounds), we can consider the range over which the marginal distribution of variance components would be sensible. For the range [-R,R], the equivalent gamma parameters:

\[ \begin{align} a &= \frac{d}{2}\\ b &= \frac{R^2d}{2*(t^d_{1-(1-q)/2})^2} \end{align} \]

where \(d\) is the degrees of freedom (for Cauchy marginal, \(d=1\)) and \(t^d_{1-(1-q)/2}\) is the \(q_{th}\) quantile of a Student t with \(d\) degrees of freedom. Hence if we considered variance likely to be in the range of [0,100], we could define a loggamma prior of:

\[ \begin{align} a &= \frac{1}{2} = 0.5\\ b &= \frac{100^2}{2*161.447} = 30.9\\ \tau &\sim{} log\Gamma(0.5,30.9) \end{align} \]

standist::visualize("gamma(0.5, 0.032)", xlim=c(-5,100))

Intercept

The default prior on the intercept is a Gaussian with mean of 0 and precision of 0 (and thus effectively a flat uniform). Alternatively, we could define priors on the intercept that are more realistic. For example, we know that the median grass yield is 80. Hence, the intercept is likely to be relatively close to this value. Note, arguably it would be more sensible to center the predictor variable (fertiliser concentration) so that the intercept would represent the yeild at the average fertiliser concentration. Furthermore, this would also be more computationally expedient. In fact, the only reason we are not doing so is so that we can directly compare the parameter estimates to the frequentist linear model.

min(fert$YIELD)
[1] 80
median(fert$YIELD)
[1] 161.5
mad(fert$YIELD)
[1] 90.4386

We will multiply the variability (90.4386) by 3 so as to estimate an envelope that so as to fully capture the full range of expected values (271.3158.

standist::visualize("normal(80, 271.32)", xlim=c(-700,1000))

We could use these values as the basis for weekly informative priors on the intercept. Note, as INLA priors are expressed in terms of precision rather than variance, an equvilent prior would be \(\sim{}~\mathcal{N}(80, 0.00001)\) (e.g. \(1/(MAD\times 3)^2\)).

Fixed effects

The priors for the fixed effects (slope) is a Gaussian (normal) distributions with mean of 0 and precision (0.001). This implies that the prior for slope has a standard deviation of approximately 31 (since \(\sigma = \frac{1}{\sqrt{\tau}}\)). As a general rule, three standard deviations envelopes most of a distribution, and thus this prior defines a distribution whose density is almost entirely within the range [-93,93]. Arguably, this is very wide for the slope given the range of the predictor variable.

In order to generate realistic informative Gaussian priors (for the purpose of constraining the posterior to a logical range) for fixed parameters, the following formulae are useful:

\[ \begin{align} \mu &= \frac{z_2\theta_1 - z_1\theta_2}{z_2-z_1}\\ \sigma &= \frac{\theta_2 - \theta_1}{z_2-z_1} \end{align} \]

where \(\theta_1\) and \(\theta_2\) are the quantiles on the response scale and \(z_1\) and \(z_2\) are the corresponding quantiles on the standard normal scale. Hence, if we considered that the slope is likely to be in the range of [-10,10], we could specify a Normal prior with mean of \(\mu=\frac{(qnorm(0.5,0,1)*0) - (qnorm(0.975,0,1)*10)}{10-0} = 0\) and a standard deviation of \(\sigma^2=\frac{10 - 0}{qnorm(0.975,0,1)-qnorm(0.5,0,1)} = 5.102\). In INLA (which defines priors in terms of precision rather than standard deviation), the associated prior would be \(\beta \sim{} \mathcal{N}(0, 0.0384)\).

standist::visualize("normal(0, 5.102)", xlim=c(-20,20))

In order to define each of the above priors, we could modify the inla call:

fert.inla1 <- inla(YIELD ~ FERTILIZER,
  data = fert,
  family = "gaussian",
  control.fixed = list(
    mean.intercept = 80,
    prec.intercept = 0.00001,
    mean = 0,
    prec = 0.0384
  ),
  control.predictor = list(compute = TRUE, link = 1),
  control.family = list(hyper = list(prec = list(
    prior = "loggamma",
    param = c(0.5, 0.032)
  ))),
  control.compute = list(
    config = TRUE,
    return.marginals.predictor = TRUE,
    dic = TRUE,
    waic = TRUE,
    cpo = TRUE
  )
)

Lets briefly compare the parameter estimates. In the following, Model.1 is the model that employed the default priors and Model.2 is the model that employed the user defined priors.

## Fixed effects
rbind(Model.1 = fert.inla$summary.fixed, Model.2 = fert.inla1$summary.fixed)
                          mean          sd 0.025quant   0.5quant 0.975quant
Model.1.(Intercept) 51.9340762 12.64763106 26.6632595 51.9339972 77.2053846
Model.1.FERTILIZER   0.8113885  0.08153394  0.6484749  0.8113891  0.9742986
Model.2.(Intercept) 52.0169173 13.48574303 25.0491429 52.0071044 79.0461379
Model.2.FERTILIZER   0.8108657  0.08695558  0.6365844  0.8109277  0.9847584
                         mode          kld
Model.1.(Intercept) 51.934013 4.155607e-06
Model.1.FERTILIZER   0.811389 4.155455e-06
Model.2.(Intercept) 52.009272 7.506730e-06
Model.2.FERTILIZER   0.810914 7.568462e-06
## Family hyperpriors (variance)
rbind(Model.1 = fert.inla$summary.hyperpar, Model.2 = fert.inla1$summary.hyperpar)
               mean          sd   0.025quant    0.5quant  0.975quant
Model.1 0.003463372 0.001543917 0.0011479418 0.003235830 0.007093335
Model.2 0.003118776 0.001463764 0.0009644277 0.002891328 0.006590397
               mode
Model.1 0.002778790
Model.2 0.002416532

It is clear that the two models yeild very similar parameter estimates.

Note, the following will not appear in this tutorial as the INLA plotting function has been really awkwardly written and is not compatible with Rmd.

plot(fert.inla1,
  plot.fixed.effects = TRUE,
  plot.lincomb = FALSE,
  plot.random.effects = FALSE,
  plot.hyperparameters = TRUE,
  plot.predictor = FALSE,
  plot.q = FALSE,
  plot.cpo = FALSE,
  plot.prior = TRUE,
  single = FALSE
)
plot_fixed_marginals(fert.inla1, priors = TRUE)

library(INLAutils)
plot_hyper_marginals(fert.inla1, CI = TRUE)

Or the effects as a single figure..

library(INLAutils)
autoplot(fert.inla1)
Warning in autoplot.inla(fert.inla1): Plot 3 selected in which, but no random
effects to plot marginals for.
Warning in autoplot.inla(fert.inla1): Plot 3 will not be plotted.

marg.fix <- fert.inla1$marginals.fixed
(nm <- names(marg.fix))
[1] "(Intercept)" "FERTILIZER" 
p <- vector("list", length(nm))
for (wch in 1:length(nm)) {
  xy <- INLA:::inla.get.prior.xy(
    section = "fixed", hyperid = nm[1],
    all.hyper = fert.inla1$all.hyper,
    range = 3 * range(marg.fix[[wch]][, "x"])
  )
  p[[wch]] <- ggplot() +
    geom_density(
      data = as.data.frame(marg.fix[[wch]]),
      aes(x = x, y = y, fill = "Posterior"), alpha = 0.2,
      stat = "Identity"
    ) +
    geom_density(
      data = as.data.frame(xy), aes(x = x, y = y, fill = "Prior"), alpha = 0.2,
      stat = "Identity"
    )
}
patchwork::wrap_plots(p)

6 MCMC sampling diagnostics

MCMC sampling behaviour

Since the purpose of the MCMC sampling is to estimate the posterior of an unknown joint likelihood, it is important that we explore a range of diagnostics designed to help identify when the resulting likelihood might not be accurate.

  • traceplots - plots of the individual draws in sequence. Traces that resemble noise suggest that all likelihood features are likely to have be traversed. Obvious steps or blocks of noise are likely to represent distinct features and could imply that there are yet other features that have not yet been traversed - necessitating additional iterations. Furthermore, each chain should be indistinguishable from the others
  • autocorrelation function - plots of the degree of correlation between pairs of draws for a range of lags (distance along the chains). High levels of correlation (after a lag of 0, which is correlating each draw with itself) suggests a lack of independence between the draws and that therefore, summaries such as mean and median will be biased estimates. Ideally, all non-zero lag correlations should be less than 0.2.
  • convergence diagnostics - there are a range of diagnostics aimed at exploring whether the multiple chains are likely to have converged upon similar posteriors
    • R hat - this metric compares between and within chain model parameter estimates, with the expectation that if the chains have converged, the between and within rank normalised estimates should be very similar (and Rhat should be close to 1). The more one chains deviates from the others, the higher the Rhat value. Values less than 1.05 are considered evidence of convergence.
    • Bulk ESS - this is a measure of the effective sample size from the whole (bulk) of the posterior and is a good measure of the sampling efficiency of draws across the entire posterior
    • Tail ESS - this is a measure of the effective sample size from the 5% and 95% quantiles (tails) of the posterior and is a good measure of the sampling efficiency of draws from the tail (areas of the posterior with least support and where samplers can get stuck).

available_mcmc()

Package Description function rstanarm brms
bayesplot Traceplot mcmc_trace plot(mod, plotfun='trace') mcmc_plot(mod, type='trace')
Density plot mcmc_dens plot(mod, plotfun='dens') mcmc_plot(mod, type='dens')
Density & Trace mcmc_combo plot(mod, plotfun='combo') mcmc_plot(mod, type='combo')
ACF mcmc_acf_bar plot(mod, plotfun='acf_bar') mcmc_plot(mod, type='acf_bar')
Rhat hist mcmc_rhat_hist plot(mod, plotfun='rhat_hist') mcmc_plot(mod, type='rhat_hist')
No. Effective mcmc_neff_hist plot(mod, plotfun='neff_hist') mcmc_plot(mod, type='neff_hist')
rstan Traceplot stan_trace stan_trace(mod) stan_trace(mod)
ACF stan_ac stan_ac(mod) stan_ac(mod)
Rhat stan_rhat stan_rhat(mod) stan_rhat(mod)
No. Effective stan_ess stan_ess(mod) stan_ess(mod)
Density plot stan_dens stan_dens(mod) stan_dens(mod)
ggmcmc Traceplot ggs_traceplot ggs_traceplot(ggs(mod)) ggs_traceplot(ggs(mod))
ACF ggs_autocorrelation ggs_autocorrelation(ggs(mod)) ggs_autocorrelation(ggs(mod))
Rhat ggs_Rhat ggs_Rhat(ggs(mod)) ggs_Rhat(ggs(mod))
No. Effective ggs_effective ggs_effective(ggs(mod)) ggs_effective(ggs(mod))
Cross correlation ggs_crosscorrelation ggs_crosscorrelation(ggs(mod)) ggs_crosscorrelation(ggs(mod))
Scale reduction ggs_grb ggs_grb(ggs(mod)) ggs_grb(ggs(mod))

In addition to the regular model diagnostics checking, for Bayesian analyses, it is also necessary to explore the MCMC sampling diagnostics to be sure that the chains are well mixed and have converged on a stable posterior.

There are a wide variety of tests that range from the big picture, overall chain characteristics to the very specific detailed tests that allow the experienced modeller to drill down to the very fine details of the chain behaviour. Furthermore, there are a multitude of packages and approaches for exploring these diagnostics.

The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.

available_mcmc()
bayesplot MCMC module:
  mcmc_acf
  mcmc_acf_bar
  mcmc_areas
  mcmc_areas_ridges
  mcmc_combo
  mcmc_dens
  mcmc_dens_chains
  mcmc_dens_overlay
  mcmc_hex
  mcmc_hist
  mcmc_hist_by_chain
  mcmc_intervals
  mcmc_neff
  mcmc_neff_hist
  mcmc_nuts_acceptance
  mcmc_nuts_divergence
  mcmc_nuts_energy
  mcmc_nuts_stepsize
  mcmc_nuts_treedepth
  mcmc_pairs
  mcmc_parcoord
  mcmc_rank_ecdf
  mcmc_rank_hist
  mcmc_rank_overlay
  mcmc_recover_hist
  mcmc_recover_intervals
  mcmc_recover_scatter
  mcmc_rhat
  mcmc_rhat_hist
  mcmc_scatter
  mcmc_trace
  mcmc_trace_highlight
  mcmc_violin

Of these, we will focus on:

  • mcmc_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
plot(fert.rstanarm3, plotfun = "mcmc_trace")

The chains appear well mixed and very similar

  • acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain
plot(fert.rstanarm3, "acf_bar")

There is no evidence of autocorrelation in the MCMC samples

  • Rhat: Rhat is a measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
plot(fert.rstanarm3, "rhat_hist")
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

  • neff (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

    If the ratios are low, tightening the priors may help.

plot(fert.rstanarm3, "neff_hist")
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios all very high.

More diagnostics
plot(fert.rstanarm3, "combo")

plot(fert.rstanarm3, "violin")

The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

  • stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
stan_trace(fert.rstanarm3)

The chains appear well mixed and very similar

  • stan_acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain
stan_ac(fert.rstanarm3)

There is no evidence of autocorrelation in the MCMC samples

  • stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
stan_rhat(fert.rstanarm3)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

  • stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

    If the ratios are low, tightening the priors may help.

stan_ess(fert.rstanarm3)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios all very high.

stan_dens(fert.rstanarm3, separate_chains = TRUE)

The ggmean package also has a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.

Of these, we will focus on:

  • ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
fert.ggs <- ggs(fert.rstanarm3)
ggs_traceplot(fert.ggs)

The chains appear well mixed and very similar

  • gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain
ggs_autocorrelation(fert.ggs)

There is no evidence of autocorrelation in the MCMC samples

  • stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
ggs_Rhat(fert.ggs)

All Rhat values are below 1.05, suggesting the chains have converged.

  • stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

    If the ratios are low, tightening the priors may help.

ggs_effective(fert.ggs)
Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
dplyr 1.1.0.
ℹ Please use `reframe()` instead.
ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
  always returns an ungrouped data frame and adjust accordingly.
ℹ The deprecated feature was likely used in the ggmcmc package.
  Please report the issue at <https://github.com/xfim/ggmcmc/issues/>.

Ratios all very high.

More diagnostics
ggs_crosscorrelation(fert.ggs)

ggs_grb(fert.ggs)

The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.

available_mcmc()
bayesplot MCMC module:
  mcmc_acf
  mcmc_acf_bar
  mcmc_areas
  mcmc_areas_ridges
  mcmc_combo
  mcmc_dens
  mcmc_dens_chains
  mcmc_dens_overlay
  mcmc_hex
  mcmc_hist
  mcmc_hist_by_chain
  mcmc_intervals
  mcmc_neff
  mcmc_neff_hist
  mcmc_nuts_acceptance
  mcmc_nuts_divergence
  mcmc_nuts_energy
  mcmc_nuts_stepsize
  mcmc_nuts_treedepth
  mcmc_pairs
  mcmc_parcoord
  mcmc_rank_ecdf
  mcmc_rank_hist
  mcmc_rank_overlay
  mcmc_recover_hist
  mcmc_recover_intervals
  mcmc_recover_scatter
  mcmc_rhat
  mcmc_rhat_hist
  mcmc_scatter
  mcmc_trace
  mcmc_trace_highlight
  mcmc_violin

Of these, we will focus on:

  • trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
fert.brm3 |> mcmc_plot(type = "combo")

fert.brm3 |> mcmc_plot(type = "trace")
No divergences to plot.

fert.brm3 |> mcmc_plot(type = "dens_overlay")

fert.brm3 |> mcmc_plot(type = "nuts_treedepth")

fert.brm3 |> mcmc_plot(type = "nuts_divergence")

fert.brm3 |> mcmc_plot(type = "nuts_stepsize")

The chains appear well mixed and very similar

  • acf_bar (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain
fert.brm3 |> mcmc_plot(type = "acf_bar")

There is no evidence of autocorrelation in the MCMC samples

  • rhat_hist: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
fert.brm3 |> mcmc_plot(type = "rhat_hist")
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

  • neff_hist (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

    If the ratios are low, tightening the priors may help.

fert.brm3 |> mcmc_plot(type = "neff_hist")
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios all very high.

More diagnostics
fert.brm3 |> mcmc_plot(type = "combo")

fert.brm3 |> mcmc_plot(type = "violin")

The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

  • stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
stan_trace(fert.brm3$fit)

fert.brm3$fit |> stan_trace()

fert.brm3$fit |> stan_trace(inc_warmup = TRUE)

The chains appear well mixed and very similar

  • stan_acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain
fert.brm3$fit |> stan_ac()

There is no evidence of autocorrelation in the MCMC samples

  • stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
fert.brm3$fit |> stan_rhat()
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

  • stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

    If the ratios are low, tightening the priors may help.

fert.brm3$fit |> stan_ess()
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios all very high.

fert.brm3$fit |> stan_dens(separate_chains = TRUE)

The ggmcmc package also has a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.

Of these, we will focus on:

  • ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
fert.ggs <- fert.brm3$fit |> ggs(inc_warmup = TRUE, burnin = TRUE) # does not seem to ignore burnin?
fert.ggs |> ggs_traceplot()

fert.ggs <- fert.brm3$fit |> ggs(inc_warmup = FALSE, burnin = TRUE) # does not seem to ignore burnin?
fert.ggs |> ggs_traceplot()

The chains appear well mixed and very similar

  • gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain
fert.ggs |> ggs_autocorrelation()

There is no evidence of autocorrelation in the MCMC samples

  • stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
fert.ggs |> ggs_Rhat()

All Rhat values are below 1.05, suggesting the chains have converged.

  • stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

    If the ratios are low, tightening the priors may help.

fert.ggs |> ggs_effective()

Ratios all very high.

More diagnostics
fert.ggs |> ggs_crosscorrelation()

fert.ggs |> ggs_grb()

The sampling diagnostics are not relevant to INLA. Instead, we can focus on CPO and PIT

fert.inla1.cpo <- data.frame(
  cpo = fert.inla1$cpo$cpo,
  pit = fert.inla1$cpo$pit,
  failure = fert.inla1$cpo$failure
) |>
  filter(failure == 0) |>
  mutate(row = 1:n())

g1 <- fert.inla1.cpo |>
  ggplot(aes(y = cpo, x = row)) +
  geom_point()
g2 <- fert.inla1.cpo |>
  ggplot(aes(x = cpo)) +
  geom_histogram()
g3 <- fert.inla1.cpo |>
  ggplot(aes(y = pit, x = row)) +
  geom_point()
g4 <- fert.inla1.cpo |>
  ggplot(aes(x = pit)) +
  geom_histogram()

(g1 + g2) / (g3 + g4)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ideally, we want to see that no CPO or PIT is very different in value from any others (not the case here) and that the histograms are relatively flat (which they kind of are here - particularly considering the small amount of data).

7 Model validation

Posterior probability checks

available_ppc()

Package Description function rstanarm brms
bayesplot Density overlay ppc_dens_overlay pp_check(mod, plotfun='dens_overlay') pp_check(mod, type='dens_overlay')
Obs vs Pred error ppc_error_scatter_avg pp_check(mod, plotfun='error_scatter_avg') pp_check(mod, type='error_scatter_avg')
Pred error vs x ppc_error_scatter_avg_vs_x pp_check(mod, x=, plotfun='error_scatter_avg_vs_x') pp_check(mod, x=, type='error_scatter_avg_vs_x')
Preds vs x ppc_intervals pp_check(mod, x=, plotfun='intervals') pp_check(mod, x=, type='intervals')
Partial plot ppc_ribbon pp_check(mod, x=, plotfun='ribbon') pp_check(mod, x=, type='ribbon')

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

available_ppc()
bayesplot PPC module:
  ppc_bars
  ppc_bars_grouped
  ppc_boxplot
  ppc_dens
  ppc_dens_overlay
  ppc_dens_overlay_grouped
  ppc_ecdf_overlay
  ppc_ecdf_overlay_grouped
  ppc_error_binned
  ppc_error_hist
  ppc_error_hist_grouped
  ppc_error_scatter
  ppc_error_scatter_avg
  ppc_error_scatter_avg_grouped
  ppc_error_scatter_avg_vs_x
  ppc_freqpoly
  ppc_freqpoly_grouped
  ppc_hist
  ppc_intervals
  ppc_intervals_grouped
  ppc_km_overlay
  ppc_km_overlay_grouped
  ppc_loo_intervals
  ppc_loo_pit_ecdf
  ppc_loo_pit_overlay
  ppc_loo_pit_qq
  ppc_loo_ribbon
  ppc_pit_ecdf
  ppc_pit_ecdf_grouped
  ppc_ribbon
  ppc_ribbon_grouped
  ppc_rootogram
  ppc_scatter
  ppc_scatter_avg
  ppc_scatter_avg_grouped
  ppc_stat
  ppc_stat_2d
  ppc_stat_freqpoly
  ppc_stat_freqpoly_grouped
  ppc_stat_grouped
  ppc_violin_grouped
  • dens_overlay: plots the density distribution of the observed data (black line) overlayed on top of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.
pp_check(fert.rstanarm3, plotfun = "dens_overlay")

The model draws appear to be consistent with the observed data.

  • error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot
pp_check(fert.rstanarm3, plotfun = "error_scatter_avg")

There is no obvious pattern in the residuals.

  • error_scatter_avg_vs_x: this is similar to a regular residual plot and as such should be interpreted as such.
pp_check(fert.rstanarm3, x = fert$FERTILIZER, plotfun = "error_scatter_avg_vs_x")

  • intervals: plots the observed data overlayed onto of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.
pp_check(fert.rstanarm3, x = fert$FERTILIZER, plotfun = "intervals")

  • ribbon: this is just an alternative way of expressing the above plot.
pp_check(fert.rstanarm3, x = fert$FERTILIZER, plotfun = "ribbon")

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

# library(shinystan)
# launch_shinystan(fert.rstanarm3)

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

  • simulated (predicted) responses associated with each observation.
  • observed values
  • fitted (predicted) responses (averaged) associated with each observation
preds <- posterior_predict(fert.rstanarm3, ndraws = 250, summary = FALSE)
fert.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = fert$YIELD,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = "gaussian"
)
plot(fert.resids)

Conclusions:

  • the simulated residuals do not suggest any issues with the residuals
  • there is no evidence of a lack of fit.

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

available_ppc()
bayesplot PPC module:
  ppc_bars
  ppc_bars_grouped
  ppc_boxplot
  ppc_dens
  ppc_dens_overlay
  ppc_dens_overlay_grouped
  ppc_ecdf_overlay
  ppc_ecdf_overlay_grouped
  ppc_error_binned
  ppc_error_hist
  ppc_error_hist_grouped
  ppc_error_scatter
  ppc_error_scatter_avg
  ppc_error_scatter_avg_grouped
  ppc_error_scatter_avg_vs_x
  ppc_freqpoly
  ppc_freqpoly_grouped
  ppc_hist
  ppc_intervals
  ppc_intervals_grouped
  ppc_km_overlay
  ppc_km_overlay_grouped
  ppc_loo_intervals
  ppc_loo_pit_ecdf
  ppc_loo_pit_overlay
  ppc_loo_pit_qq
  ppc_loo_ribbon
  ppc_pit_ecdf
  ppc_pit_ecdf_grouped
  ppc_ribbon
  ppc_ribbon_grouped
  ppc_rootogram
  ppc_scatter
  ppc_scatter_avg
  ppc_scatter_avg_grouped
  ppc_stat
  ppc_stat_2d
  ppc_stat_freqpoly
  ppc_stat_freqpoly_grouped
  ppc_stat_grouped
  ppc_violin_grouped
  • dens_overlay: plots the density distribution of the observed data (black line) overlayed on top of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.
fert.brm3 |> pp_check(type = "dens_overlay", ndraws = 100)

The model draws appear to be consistent with the observed data.

  • error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot
fert.brm3 |> pp_check(type = "error_scatter_avg")
Using all posterior draws for ppc type 'error_scatter_avg' by default.

There is no obvious pattern in the residuals.

  • error_scatter_avg_vs_x: this is similar to a regular residual plot and as such should be interpreted as such.
fert.brm3 |> pp_check(x = "FERTILIZER", type = "error_scatter_avg_vs_x")
Using all posterior draws for ppc type 'error_scatter_avg_vs_x' by default.

  • intervals: plots the observed data overlayed on top of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.
fert.brm3 |> pp_check(x = "FERTILIZER", type = "intervals")
Using all posterior draws for ppc type 'intervals' by default.

  • ribbon: this is just an alternative way of expressing the above plot.
pp_check(fert.brm3, x = "FERTILIZER", type = "ribbon")
Using all posterior draws for ppc type 'ribbon' by default.

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

# library(shinystan)
# launch_shinystan(fert.brm3)

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

  • simulated (predicted) responses associated with each observation.
  • observed values
  • fitted (predicted) responses (averaged) associated with each observation
preds <- fert.brm3 |> posterior_predict(ndraws = 250, summary = FALSE)
fert.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = fert$YIELD,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = FALSE
)
fert.resids |> plot()

fert.resids <- make_brms_dharma_res(fert.brm3, integerResponse = FALSE)
wrap_elements(~ testUniformity(fert.resids)) +
  wrap_elements(~ plotResiduals(fert.resids, form = factor(rep(1, nrow(fert))))) +
  wrap_elements(~ plotResiduals(fert.resids)) +
  wrap_elements(~ testDispersion(fert.resids))

The integerResponse argument indicates that noise will be added to the residuals such that they will become integers. This is important in for Poisson, Negative Binomial and Binomial models.

Conclusions:

  • the simulated residuals do not suggest any issues with the residuals
  • there is no evidence of a lack of fit.

INLA includes two measures of out-of-sample performace (measures of fit).

  • Conditional Predictive Ordinate (CPO) is the density of the observed value of \(y_i\) within the out-of-sample (\(y−i\)) posterior predictive distribution. A small CPO value associated with an observation suggests that this observation is unlikely (or surprising) in light of the model, priors and other data in the model. In addition, the sum of the CPO values (or alternatively, the negative of the mean natural logarithm of the CPO values) is a measure of fit.

    fert.inla1$cpo$cpo
     [1] 0.012299827 0.013183921 0.006748973 0.008780066 0.018831751 0.019056236
     [7] 0.015164579 0.002741200 0.005978481 0.015229534
    sum(fert.inla1$cpo$cop)
    [1] 0
    -mean(log(fert.inla1$cpo$cpo))
    [1] 4.57854

    In this case there are no CPO values that are considerably smaller (order of magnitude smaller) than the others - so with respect to the model, none of the observed values would be considered ‘surprising’. Various assumptions are made behind INLA’s computations. If any of these assumptions fail for an observation, it is flagged within (in this case) cpo$failure as a non-zero. When this is the case, it is prudent to recalculate the CPO values after removing the observations associated with failure flags not equal to zero and re-fitting the model. This can be performed using the inla.cpo() function.

    fert.inla1$cpo$failure
     [1] 0 0 0 0 0 0 0 0 0 0

    In this case, there are no non-zero CPO values.

  • Probability Integral Transforms (PIT) provides a version of CPO that is calibrated to the level of the Gaussian field so that it is clearer whether or not any of the values are ‘small’ (all values must be between 0 and 1). A histogram of PIT values that does not look approximately uniform (flat) indicates a lack of model fit.

    Unfortunately, the following will not work inside of an Rmarkdown (but will work fine on the console)

In this case the PIT values do not appear to deviate from a uniform distribution and thus do not indicate a lack of fit.

fert.inla1.cpo <- data.frame(
  cpo = fert.inla1$cpo$cpo,
  pit = fert.inla1$cpo$pit,
  failure = fert.inla1$cpo$failure
) |>
  filter(failure == 0) |>
  mutate(row = 1:n())

g1 <- fert.inla1.cpo |>
  ggplot(aes(y = cpo, x = row)) +
  geom_point()
g2 <- fert.inla1.cpo |>
  ggplot(aes(x = cpo)) +
  geom_histogram()
g3 <- fert.inla1.cpo |>
  ggplot(aes(y = pit, x = row)) +
  geom_point()
g4 <- fert.inla1.cpo |>
  ggplot(aes(x = pit)) +
  geom_histogram()

(g1 + g2) / (g3 + g4)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ideally, we want to see that no CPO or PIT is very different in value from any others (not the case here) and that the histograms are relatively flat (which they kind of are here - particularly considering the small amount of data).

In the following, the top left is the distribution of residuals and top right is the observed vs predicted values.

ggplot_inla_residuals(fert.inla1, observed = fert$YIELD)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Conclusions-

  • the distribution of residuals is nearly uniform (keep in mind that the sample size is very small in this example).
  • the fitted values approximate the observed values very well

The following figure represents the equivalent of a (standardised) residual plot.

ggplot_inla_residuals2(fert.inla1, observed = fert$YIELD)
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Unfortunately, DHARMa does not directly support INLA. Therefore, in order to make use of this framework, we need to do some of the heavy lifting. Specifically, we need to generate the posterior predictions associated with the observed data.

Whilst we can draw random samples from the posterior with the inla.posterior.sample() function, this will only account for uncertainty in the fixed parameter estimates. That is, it will not incorporate the overall uncertainty in the Gaussian distribution (\(\sigma\)). Therefore, we need to add this back on ourself (keep in mind that this will be in units of precision and thus needs to be back-scaled into units of standard deviation). Furthermore, we will need to use the estimates of standard deviation to estimate the random noise to add onto the predictions. This is done by resampling the family (Gaussian in this case).

I acknowledge that this (resampling the SD) appears to be a bit of a hack and seems to go against the INLA principles. Nevertheless, we are only doing this this for the purposes of checking residual diagnostics.

## draw 250 samples from the posterior
draws <- inla.posterior.sample(n = 250, result = fert.inla)
## extract the latent predictions for the observed data
preds <- t(sapply(draws, function(x) x$latent[1:nrow(fert)]))
## extract the first (family precision) hyperprior
preds.hyper <- sapply(draws, function(x) 1 / sqrt(x$hyperpar[[1]]))
## use this to generate gaussian noise
preds.hyper <- rnorm(length(preds.hyper), mean = 0, sd = preds.hyper)
## add the noise to each prediction
preds <- sweep(preds, 1, preds.hyper, FUN = "+")
## generate the DHARMa residuals
fert.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = fert$YIELD,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = FALSE
)
fert.resids |> plot()

preds <- posterior_predict.inla(fert.inla, newdata = fert)
fert.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = fert$YIELD,
  fittedPredictedResponse = apply(preds, 2, mean),
  integerResponse = FALSE
)
wrap_elements(~ testUniformity(fert.resids)) +
  wrap_elements(~ plotResiduals(fert.resids, form = factor(rep(1, nrow(fert))))) +
  wrap_elements(~ plotResiduals(fert.resids)) +
  wrap_elements(~ testDispersion(fert.resids))

Conclusions:

  • the simulated residuals do not suggest any issues with the residuals
  • there is no evidence of a lack of fit.

8 Partial effects plots

fert.rstanarm3 |>
  ggpredict() |>
  plot(show_data = TRUE)
Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

fert.rstanarm3 |>
  ggemmeans(~FERTILIZER) |>
  plot(show_data = TRUE)
Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

fert.rstanarm3 |>
  epred_draws(newdata = fert) |>
  median_hdci() |>
  ggplot(aes(x = FERTILIZER, y = .epred)) +
  geom_ribbon(aes(ymin = .lower, ymax = .upper), fill = "blue", alpha = 0.3) +
  geom_line()

fert.brm3 |> conditional_effects()

fert.brm3 |>
  conditional_effects() |>
  plot(points = TRUE)

fert.brm3 |>
  conditional_effects(spaghetti = TRUE, ndraws = 200) |>
  plot(points = TRUE)

fert.brm3 |>
  ggpredict() |>
  plot(show_data = TRUE)
Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

fert.brm3 |>
  ggemmeans(~FERTILIZER) |>
  plot(show_data = TRUE)
Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

fert.brm3 |>
  epred_draws(newdata = fert) |>
  median_hdci() |>
  ggplot(aes(x = FERTILIZER, y = .epred)) +
  geom_ribbon(aes(ymin = .lower, ymax = .upper), fill = "blue", alpha = 0.3) +
  geom_line()

plot(fert.inla1,
  plot.fixed.effects = FALSE,
  plot.lincomb = FALSE,
  plot.random.effects = FALSE,
  plot.hyperparameters = FALSE,
  plot.predictor = TRUE,
  plot.q = FALSE,
  plot.cpo = FALSE,
  plot.prior = FALSE
)
newdata <- fert |> tidyr::expand(FERTILIZER = modelr::seq_range(FERTILIZER, n = 100))
Xmat <- model.matrix(~FERTILIZER, newdata)

nms <- colnames(fert.inla1$model.matrix)
n <- sapply(nms, function(x) 0, simplify = FALSE)
draws <- inla.posterior.sample(n = 250, result = fert.inla1, selection = n)
coefs <- t(sapply(draws, function(x) x$latent))
Fit <- coefs %*% t(Xmat) |>
  as.data.frame() |>
  mutate(Rep = 1:n()) |>
  pivot_longer(cols = -Rep) |>
  group_by(name) |>
  median_hdci(value) |>
  ungroup() |>
  mutate(name = as.integer(as.character(name))) |>
  arrange(name)
newdata <- newdata |>
  bind_cols(Fit)

newdata |>
  ggplot(aes(x = FERTILIZER)) +
  geom_ribbon(aes(ymin = .lower, ymax = .upper), fill = "orange", color = NA, alpha = 0.3) +
  geom_line(aes(y = value), color = "orange") +
  geom_point(data = fert, aes(y = YIELD)) +
  theme_classic()

9 Model investigation

Rather than simply return point estimates of each of the model parameters, Bayesian analyses capture the full posterior of each parameter. These are typically stored within the list structure of the output object.

As with most statistical routines, the overloaded summary() function provides an overall summary of the model parameters. Typically, the summaries will include the means / medians along with credibility intervals and perhaps convergence diagnostics (such as R hat). However, more thorough investigation and analysis of the parameter posteriors requires access to the full posteriors.

There is currently a plethora of functions for extracting the full posteriors from models. In part, this is a reflection of a rapidly evolving space with numerous packages providing near equivalent functionality (it should also be noted, that over time, many of the functions have been deprecated due to inconsistencies in their names). Broadly speaking, the functions focus on draws from the posterior of either the parameters (intercept, slope, standard deviation etc), linear predictor, expected values or predicted values. The distinction between the latter three are highlighted in the following table.

Property Description
linear predictors values predicted on the link scale
expected values predictions (on response scale) without residual error (predicting expected mean outcome(s))
predicted values predictions (on response scale) that incorporate residual error
fitted values predictions on the response scale

The following table lists the various ways of extracting the full posteriors of the model parameters parameters, expected values and predicted values. The crossed out items are now deprecated and function with a namespace of __ mean that the functionality is provided via a range of packages.

Function Values Description
__::as.matrix() Parameters Returns \(n\times p\) matrix
__::as.data.frame() Parameters Returns \(n\times p\) data.frame
__::as_tibble() Parameters Returns \(n\times p\) tibble
posterior::as_draws_df() Parameters Returns \(n\times p\) data.frame with additional info about chain, interaction and draw
brms::posterior_samples() Parameters Returns \(n\times p\) data.frame
tidybayes::tidy_draws() Parameters Returns \(n\times p\) tibble with addition info about the chain, iteration and draw
rstan::extract() Parameters Returns a \(p\) length list of \(n\) length vectors
tidybayes::spread_draws() Parameters Returns \(n\times r\) tibble with additional info about chain, interaction and draw
tidybayes::gather_draws() Parameters Returns a gathered spread_draws tibble with additional info about chain, interaction and draw
rstanarm::posterior_linpred() Linear predictors Returns \(n\times N\) tibble on the link scale
brms::posterior_linpred() Linear predictors Returns \(n\times N\) tibble on the link scale
tidybayes::linpred_draws() Linear predictors Returns tibble with $nN rows and .linpred on the link scale additional info about chain, interaction and draw
rstanarm::posterior_epred() Expected values Returns \(n\times N\) tibble on the response scale
brms::posterior_epred() Expected values Returns \(n\times N\) tibble on the response scale
tidybayes::epred_draws() Expected values Returns tibble with $nN rows and .epred on the response scale additional info about chain, interaction and draw
rstanarm::posterior_predict() Expected values Returns \(n\times N\) tibble of predictions (including residuals) on the response scale
brms::posterior_predict() Expected values Returns \(n\times N\) tibble of predictions (including residuals) on the response scale
tidybayes::predicted_draws() Expected values Returns tibble with $nN rows and .prediction (including residuals) on the response scale additional info about chain, interaction and draw

where \(n\) is the number of MCMC samples and \(p\) is the number of parameters to estimate, \(N\) is the number of newdata rows and \(r\) is the number of requested parameters. For the tidybayes versions in the table above, the function expects a model to be the first parameter (and a dataframe to be the second). There are also add_ versions which expect a dataframe to be the first argument and the model to be the second. These alternatives facilitate pipings with different starting objects.

Function Description
median_qi Median and quantiles of specific columns
median_hdi Median and Highest Probability Density Interval of specific columns
median_hdci Median and continuous Highest Probability Density Interval of specific columns
tidyMCMC Median/mean and quantiles/hpd of all columns

It is important to remember that when working with Bayesian models, everything is a distribution. Whilst point estimates (such as a mean) of the parameters can be calculated from these distributions, we start off with a large number of estimates per parameter.

rstanarm captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

fert.rstanarm3 |> summary()

Model Info:
 function:     stan_glm
 family:       gaussian [identity]
 formula:      YIELD ~ FERTILIZER
 algorithm:    sampling
 sample:       2400 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 10
 predictors:   2

Estimates:
              mean   sd   10%   50%   90%
(Intercept) 52.3   15.7 33.4  52.0  71.9 
FERTILIZER   0.8    0.1  0.7   0.8   0.9 
sigma       22.6    7.0 15.5  21.3  31.5 

Fit Diagnostics:
           mean   sd    10%   50%   90%
mean_PPD 163.5   10.5 150.9 163.4 176.2

The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).

MCMC diagnostics
              mcse Rhat n_eff
(Intercept)   0.3  1.0  2360 
FERTILIZER    0.0  1.0  2556 
sigma         0.2  1.0  2147 
mean_PPD      0.2  1.0  2212 
log-posterior 0.0  1.0  2150 

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Conclusions:

  • in the Model Info, we are informed that the total MCMC posterior sample size is 2400 and that there were 10 raw observations.
  • the estimated intercept (expected yield when fertilizer concentration is 0) is 52.28. This is the mean of the posterior distribution for this parameter.
  • the estimated slope (rate at which yield changes per 1 unit change in fertilizer concentration), is 0.81 (mean) or 0.69 (median) with a standard deviation of 0. The 90% credibility intervals indicate that we are 90% confident that the slope is between 0.81 and 0.81 - e.g. there is a significant positive trend.
  • sigma is estimated to be 22.64
  • Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.
tidyMCMC(fert.rstanarm3$stanfit,
  estimate.method = "median", conf.int = TRUE,
  conf.method = "HPDinterval", rhat = TRUE, ess = TRUE
)
# A tibble: 5 × 7
  term          estimate std.error conf.low conf.high  rhat   ess
  <chr>            <dbl>     <dbl>    <dbl>     <dbl> <dbl> <int>
1 (Intercept)     52.3      15.7     22.0       84.3  1.00   2360
2 FERTILIZER       0.808     0.101    0.616      1.02 1.00   2556
3 sigma           22.6       7.02    11.5       36.2  1.00   2147
4 mean_PPD       163.       10.5    143.       184.   0.999  2212
5 log-posterior  -52.1       1.43   -54.8      -50.3  1.000  2150

Conclusions:

  • the estimated intercept (expected yield when fertilizer concentration is 0) is 52.28. This is the median of the posterior distribution for this parameter.
  • the estimated slope (rate at which yield changes per 1 unit change in fertilizer concentration), is 0.81 (median) with a standard error of 0.1. The 95% credibility intervals indicate that we are 95% confident that the slope is between 0.62 and 1.02 - e.g. there is a significant positive trend.
  • sigma is estimated to be 22.64
  • Rhat and number of effective samples for each parameter are also provided and all look good.
fert.rstanarm3$stanfit |>
  summarise_draws(
    median,
    HDInterval::hdi,
    rhat, length, ess_bulk, ess_tail
  )
# A tibble: 5 × 8
  variable       median   lower  upper  rhat length ess_bulk ess_tail
  <chr>           <dbl>   <dbl>  <dbl> <dbl>  <dbl>    <dbl>    <dbl>
1 (Intercept)    52.0    22.0    84.3  1.00    2400    2379.    2129.
2 FERTILIZER      0.809   0.616   1.02 1.00    2400    2585.    2500.
3 sigma          21.3    11.5    36.2  1.00    2400    2228.    2181.
4 mean_PPD      163.    143.    184.   1.000   2400    2223.    2106.
5 log-posterior -51.7   -54.8   -50.3  1.00    2400    2155.    2110.

We can also alter the CI level.

fert.rstanarm3$stanfit |>
  summarise_draws(
    median,
    ~ HDInterval::hdi(.x, credMass = 0.9),
    rhat, length, ess_bulk, ess_tail
  )
# A tibble: 5 × 8
  variable       median   lower   upper  rhat length ess_bulk ess_tail
  <chr>           <dbl>   <dbl>   <dbl> <dbl>  <dbl>    <dbl>    <dbl>
1 (Intercept)    52.0    26.7    76.7   1.00    2400    2379.    2129.
2 FERTILIZER      0.809   0.642   0.966 1.00    2400    2585.    2500.
3 sigma          21.3    12.7    32.6   1.00    2400    2228.    2181.
4 mean_PPD      163.    146.    180.    1.000   2400    2223.    2106.
5 log-posterior -51.7   -53.9   -50.3   1.00    2400    2155.    2110.
fert.rstanarm3$stanfit |> as_draws_df()
# A draws_df: 800 iterations, 3 chains, and 5 variables
   (Intercept) FERTILIZER sigma mean_PPD log-posterior
1           38       0.86    18      149           -51
2           49       0.82    17      162           -50
3           31       0.93    26      158           -52
4           52       0.77    32      151           -53
5           42       0.83    13      158           -53
6           29       0.92    36      143           -54
7           84       0.65    13      179           -58
8           32       0.86    29      149           -53
9           97       0.65    31      178           -56
10          83       0.65    18      170           -53
# ... with 2390 more draws
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}
fert.rstanarm3$stanfit |>
  as_draws_df() |>
  summarise_draws(
    median,
    ~ HDInterval::hdi(.x),
    rhat,
    ess_bulk, ess_tail
  )
# A tibble: 5 × 7
  variable       median   lower  upper  rhat ess_bulk ess_tail
  <chr>           <dbl>   <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 (Intercept)    52.0    22.0    84.3  1.00     2379.    2129.
2 FERTILIZER      0.809   0.616   1.02 1.00     2585.    2500.
3 sigma          21.3    11.5    36.2  1.00     2228.    2181.
4 mean_PPD      163.    143.    184.   1.000    2223.    2106.
5 log-posterior -51.7   -54.8   -50.3  1.00     2155.    2110.

This is purely a graphical depiction on the posteriors.

fert.rstanarm3 |> tidy_draws()
# A tibble: 2,400 × 12
   .chain .iteration .draw `(Intercept)` FERTILIZER sigma accept_stat__
    <int>      <int> <int>         <dbl>      <dbl> <dbl>         <dbl>
 1      1          1     1          38.3      0.865  18.4         1    
 2      1          2     2          49.2      0.818  17.3         0.833
 3      1          3     3          31.0      0.929  25.7         0.999
 4      1          4     4          52.5      0.772  31.7         0.994
 5      1          5     5          42.3      0.830  13.2         0.941
 6      1          6     6          29.1      0.918  35.6         0.605
 7      1          7     7          83.7      0.645  13.0         0.816
 8      1          8     8          32.5      0.865  28.8         0.973
 9      1          9     9          97.2      0.650  31.0         0.980
10      1         10    10          82.9      0.651  18.5         0.915
# ℹ 2,390 more rows
# ℹ 5 more variables: stepsize__ <dbl>, treedepth__ <dbl>, n_leapfrog__ <dbl>,
#   divergent__ <dbl>, energy__ <dbl>
fert.draw <- fert.rstanarm3 |> gather_draws(`(Intercept)`, FERTILIZER, sigma)
## OR via regex
fert.draw <- fert.rstanarm3 |> gather_draws(`.Intercept.*|FERT.*|sigma`, regex = TRUE)
fert.draw
# A tibble: 7,200 × 5
# Groups:   .variable [3]
   .chain .iteration .draw .variable   .value
    <int>      <int> <int> <chr>        <dbl>
 1      1          1     1 (Intercept)   38.3
 2      1          2     2 (Intercept)   49.2
 3      1          3     3 (Intercept)   31.0
 4      1          4     4 (Intercept)   52.5
 5      1          5     5 (Intercept)   42.3
 6      1          6     6 (Intercept)   29.1
 7      1          7     7 (Intercept)   83.7
 8      1          8     8 (Intercept)   32.5
 9      1          9     9 (Intercept)   97.2
10      1         10    10 (Intercept)   82.9
# ℹ 7,190 more rows

We can then summarise this

fert.draw |> median_hdci()
# A tibble: 3 × 7
  .variable   .value .lower .upper .width .point .interval
  <chr>        <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 (Intercept) 52.0   22.0    84.3    0.95 median hdci     
2 FERTILIZER   0.809  0.619   1.02   0.95 median hdci     
3 sigma       21.3   11.7    36.7    0.95 median hdci     
fert.rstanarm3 |>
  gather_draws(`(Intercept)`, FERTILIZER, sigma) |>
  ggplot() +
  stat_halfeye(aes(x = .value, y = .variable)) +
  facet_wrap(~.variable, scales = "free")

Conclusions:

  • the estimated intercept (expected yield when fertilizer concentration is 0) is 52.04. This is the median of the posterior distribution for this parameter.
  • the estimated slope (rate at which yield changes per 1 unit change in fertilizer concentration), is 0.81 (mean) or (median). The 95% credibility intervals indicate that we are 95% confident that the slope is between 0.62 and 1.02 - e.g. there is a significant positive trend.
  • sigma is estimated to be 21.26
fert.rstanarm3 |> plot(plotfun = "mcmc_intervals")

fert.rstanarm3 |> spread_draws(`(Intercept)`, FERTILIZER, sigma)
# A tibble: 2,400 × 6
   .chain .iteration .draw `(Intercept)` FERTILIZER sigma
    <int>      <int> <int>         <dbl>      <dbl> <dbl>
 1      1          1     1          38.3      0.865  18.4
 2      1          2     2          49.2      0.818  17.3
 3      1          3     3          31.0      0.929  25.7
 4      1          4     4          52.5      0.772  31.7
 5      1          5     5          42.3      0.830  13.2
 6      1          6     6          29.1      0.918  35.6
 7      1          7     7          83.7      0.645  13.0
 8      1          8     8          32.5      0.865  28.8
 9      1          9     9          97.2      0.650  31.0
10      1         10    10          82.9      0.651  18.5
# ℹ 2,390 more rows
# OR via regex
fert.rstanarm3 |> spread_draws(`.Intercept.*|FERT.*|sigma`, regex = TRUE)
# A tibble: 2,400 × 6
   .chain .iteration .draw `(Intercept)` FERTILIZER sigma
    <int>      <int> <int>         <dbl>      <dbl> <dbl>
 1      1          1     1          38.3      0.865  18.4
 2      1          2     2          49.2      0.818  17.3
 3      1          3     3          31.0      0.929  25.7
 4      1          4     4          52.5      0.772  31.7
 5      1          5     5          42.3      0.830  13.2
 6      1          6     6          29.1      0.918  35.6
 7      1          7     7          83.7      0.645  13.0
 8      1          8     8          32.5      0.865  28.8
 9      1          9     9          97.2      0.650  31.0
10      1         10    10          82.9      0.651  18.5
# ℹ 2,390 more rows

Note, this is not deprecated

fert.rstanarm3 |>
  posterior_samples() |>
  as_tibble()
Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.
# A tibble: 2,400 × 3
   `(Intercept)` FERTILIZER sigma
           <dbl>      <dbl> <dbl>
 1          38.3      0.865  18.4
 2          49.2      0.818  17.3
 3          31.0      0.929  25.7
 4          52.5      0.772  31.7
 5          42.3      0.830  13.2
 6          29.1      0.918  35.6
 7          83.7      0.645  13.0
 8          32.5      0.865  28.8
 9          97.2      0.650  31.0
10          82.9      0.651  18.5
# ℹ 2,390 more rows
fert.rstanarm3 |>
  bayes_R2() |>
  median_hdci()
         y      ymin      ymax .width .point .interval
1 0.892953 0.7197469 0.9704468   0.95 median      hdci

It is possible to obtain an empirical p-value (for those that cannot live without p-values). This is essentially compares (for any posterior) that the mean is zero compared to a multivariate distribution with elliptical contours.

mcmcpvalue <- function(samp) {
  ## elementary version that creates an empirical p-value for the
  ## hypothesis that the columns of samp have mean zero versus a
  ## general multivariate distribution with elliptical contours.

  ## differences from the mean standardized by the observed
  ## variance-covariance factor

  ## Note, I put in the bit for single terms
  if (length(dim(samp)) == 0) {
    std <- backsolve(chol(var(samp)), cbind(0, t(samp)) - mean(samp), transpose = TRUE)
    sqdist <- colSums(std * std)
    sum(sqdist[-1] > sqdist[1]) / length(samp)
  } else {
    std <- backsolve(chol(var(samp)), cbind(0, t(samp)) - colMeans(samp), transpose = TRUE)
    sqdist <- colSums(std * std)
    sum(sqdist[-1] > sqdist[1]) / nrow(samp)
  }
}

mcmcpvalue(as.matrix(fert.rstanarm3)[, c("FERTILIZER")])
[1] 0

brms captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

fert.brm3 |> summary()
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: YIELD ~ scale(FERTILIZER, scale = FALSE) 
   Data: fert (Number of observations: 10) 
  Draws: 3 chains, each with iter = 5000; warmup = 1000; thin = 5;
         total post-warmup draws = 2400

Regression Coefficients:
                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept                     163.41      7.35   148.73   178.64 1.00     2540
scaleFERTILIZERscaleEQFALSE     0.81      0.10     0.60     1.01 1.00     2323
                            Tail_ESS
Intercept                       2325
scaleFERTILIZERscaleEQFALSE     2109

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    22.55      7.02    13.41    39.22 1.00     2333     2269

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Conclusions:

  • in the initial block of information, we are reminded on the formula as well as chain fitting characteristics. We are also informed that the total number of post-warmup MCMC samples is is 2400.
  • the estimated intercept (expected yield when fertilizer concentration is 0) is 163.41. This is the mean of the posterior distribution for this parameter.
  • the estimated slope (rate at which yield changes per 1 unit change in fertilizer concentration), is 0.81 (mean) with a standard error of 0.1. The 95% credibility intervals indicate that we are 95% confident that the slope is between 0.6 and 1.01 - e.g. there is a significant positive trend.
  • sigma is estimated to be 22.55
  • Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.
fert.brm3$fit |>
  tidyMCMC(
    estimate.method = "median",
    conf.int = TRUE,
    conf.method = "HPDinterval",
    rhat = TRUE, ess = TRUE
  )
# A tibble: 8 × 7
  term                         estimate std.error conf.low conf.high  rhat   ess
  <chr>                           <dbl>     <dbl>    <dbl>     <dbl> <dbl> <int>
1 b_Intercept                  163.        7.35   150.        179.   1.00   2456
2 b_scaleFERTILIZERscaleEQFAL…   0.805     0.102    0.599       1.00 1.000  2314
3 sigma                         22.5       7.02    11.9        35.6  1.00   2215
4 Intercept                    163.        7.35   150.        179.   1.00   2456
5 prior_Intercept              160.       87.8    -12.6       331.   1.00   2188
6 prior_b                       -0.0235    2.02    -3.91        3.96 1.00   2352
7 prior_sigma                   97.7     114.       0.0198    283.   1.000  2431
8 lprior                       -12.0       0.0407 -12.1       -11.9  1.000  2132

Conclusions:

  • the estimated intercept (expected yield when fertilizer concentration is 0) is 163.41. This is the median of the posterior distribution for this parameter.
  • the estimated slope (rate at which yield changes per 1 unit change in fertilizer concentration), is 0.81 (median) with a standard error of 0.1. The 95% credibility intervals indicate that we are 95% confident that the slope is between 0.6 and 1 - e.g. there is a significant positive trend.
  • sigma is estimated to be 22.55
  • Rhat and number of effective samples for each parameter are also provided and all look good.
fert.brm3 |>
  summarise_draws(
    median,
    HDInterval::hdi,
    rhat, length, ess_bulk, ess_tail
  )
# A tibble: 9 × 8
  variable                 median    lower  upper  rhat length ess_bulk ess_tail
  <chr>                     <dbl>    <dbl>  <dbl> <dbl>  <dbl>    <dbl>    <dbl>
1 b_Intercept            163.     150.     179.   1.00    2400    2540.    2325.
2 b_scaleFERTILIZERscal…   0.806    0.599    1.00 1.00    2400    2323.    2109.
3 sigma                   21.3     11.9     35.6  1.00    2400    2333.    2269.
4 Intercept              163.     150.     179.   1.00    2400    2540.    2325.
5 prior_Intercept        161.     -12.6    331.   1.00    2400    2232.    2293.
6 prior_b                 -0.0607  -3.91     3.96 1.00    2400    2356.    2378.
7 prior_sigma             68.8      0.0198 283.   0.999   2400    2394.    2369.
8 lprior                 -12.0    -12.1    -11.9  1.00    2400    2288.    2096.
9 lp__                   -53.0    -56.1    -51.6  1.00    2400    2227.    2012.

We can also alter the CI level.

fert.brm3 |>
  summarise_draws(
    median,
    ~ HDInterval::hdi(.x, credMass = 0.9),
    rhat, length, ess_bulk, ess_tail
  )
# A tibble: 9 × 8
  variable                 median   lower   upper  rhat length ess_bulk ess_tail
  <chr>                     <dbl>   <dbl>   <dbl> <dbl>  <dbl>    <dbl>    <dbl>
1 b_Intercept            163.     152.    175.    1.00    2400    2540.    2325.
2 b_scaleFERTILIZERscal…   0.806    0.628   0.964 1.00    2400    2323.    2109.
3 sigma                   21.3     12.2    31.4   1.00    2400    2333.    2269.
4 Intercept              163.     152.    175.    1.00    2400    2540.    2325.
5 prior_Intercept        161.      15.2   299.    1.00    2400    2232.    2293.
6 prior_b                 -0.0607  -3.29    3.24  1.00    2400    2356.    2378.
7 prior_sigma             68.8      0.263 212.    0.999   2400    2394.    2369.
8 lprior                 -12.0    -12.0   -11.9   1.00    2400    2288.    2096.
9 lp__                   -53.0    -55.3   -51.6   1.00    2400    2227.    2012.

To narrow down to just the parameters of interest, see the code under the tidy_draws tab.

fert.brm3 |> as_draws_df()
# A draws_df: 800 iterations, 3 chains, and 9 variables
   b_Intercept b_scaleFERTILIZERscaleEQFALSE sigma Intercept prior_Intercept
1          166                          0.89    16       166             213
2          167                          0.95    23       167             135
3          161                          0.80    15       161             293
4          162                          0.83    16       162             241
5          163                          0.78    19       163             152
6          161                          0.78    15       161             133
7          167                          0.83    21       167              77
8          163                          0.87    27       163             144
9          160                          0.83    23       160             107
10         140                          0.58    26       140             220
   prior_b prior_sigma lprior
1   -2.090         266    -12
2    1.385          62    -12
3    1.674          95    -12
4   -4.032          47    -12
5    0.221          65    -12
6    0.034           6    -12
7   -2.712          13    -12
8    0.805          76    -12
9   -1.147         136    -12
10   1.702         371    -12
# ... with 2390 more draws, and 1 more variables
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}

Return the draws (samples) for each parameter in wide format

fert.brm3 |> tidy_draws()
# A tibble: 2,400 × 18
   .chain .iteration .draw b_Intercept b_scaleFERTILIZERscaleE…¹ sigma Intercept
    <int>      <int> <int>       <dbl>                     <dbl> <dbl>     <dbl>
 1      1          1     1        166.                     0.890  16.2      166.
 2      1          2     2        167.                     0.954  23.2      167.
 3      1          3     3        161.                     0.796  14.8      161.
 4      1          4     4        162.                     0.826  16.1      162.
 5      1          5     5        163.                     0.777  18.6      163.
 6      1          6     6        161.                     0.778  14.7      161.
 7      1          7     7        167.                     0.833  20.6      167.
 8      1          8     8        163.                     0.871  26.8      163.
 9      1          9     9        160.                     0.826  23.2      160.
10      1         10    10        140.                     0.584  26.4      140.
# ℹ 2,390 more rows
# ℹ abbreviated name: ¹​b_scaleFERTILIZERscaleEQFALSE
# ℹ 11 more variables: prior_Intercept <dbl>, prior_b <dbl>, prior_sigma <dbl>,
#   lprior <dbl>, lp__ <dbl>, accept_stat__ <dbl>, stepsize__ <dbl>,
#   treedepth__ <dbl>, n_leapfrog__ <dbl>, divergent__ <dbl>, energy__ <dbl>
fert.brm3 |> get_variables()
 [1] "b_Intercept"                   "b_scaleFERTILIZERscaleEQFALSE"
 [3] "sigma"                         "Intercept"                    
 [5] "prior_Intercept"               "prior_b"                      
 [7] "prior_sigma"                   "lprior"                       
 [9] "lp__"                          "accept_stat__"                
[11] "stepsize__"                    "treedepth__"                  
[13] "n_leapfrog__"                  "divergent__"                  
[15] "energy__"                     
fert.brm3$fit |>
  tidy_draws() |>
  dplyr::select(matches("^b_.*|^sigma")) |>
  summarise_draws(
    median,
    HDInterval::hdi,
    rhat, length, ess_bulk, ess_tail
  )
# A tibble: 3 × 8
  variable                  median   lower  upper  rhat length ess_bulk ess_tail
  <chr>                      <dbl>   <dbl>  <dbl> <dbl>  <dbl>    <dbl>    <dbl>
1 b_Intercept              163.    150.    179.   1.000   2400    2521.    2315.
2 b_scaleFERTILIZERscaleE…   0.806   0.599   1.00 1.000   2400    2311.    2106.
3 sigma                     21.3    11.9    35.6  1.000   2400    2326.    2253.
## if we want to express the slope in terms of the original scale of the
## fertilizer concentration, we need to back scale the slope
## fert.brm3$data
## attr(fert.brm3$data[, 2], "scaled:scale")
# fert_sigma <- attr(fert.brm3$data[, "scale(FERTILIZER)"], "scaled:scale")
fert.brm3$fit |>
  tidy_draws() |>
  dplyr::select(matches("^b_.*|^sigma")) |>
  ## mutate(FERTILIZER =  b_scaleFERTILIZER / fert_sigma)  |>
  summarise_draws(
    median,
    HDInterval::hdi,
    rhat, length, ess_bulk, ess_tail
  )
# A tibble: 3 × 8
  variable                  median   lower  upper  rhat length ess_bulk ess_tail
  <chr>                      <dbl>   <dbl>  <dbl> <dbl>  <dbl>    <dbl>    <dbl>
1 b_Intercept              163.    150.    179.   1.000   2400    2521.    2315.
2 b_scaleFERTILIZERscaleE…   0.806   0.599   1.00 1.000   2400    2311.    2106.
3 sigma                     21.3    11.9    35.6  1.000   2400    2326.    2253.
## for just the slope...
## fert.brm3 |> emtrends(~ 0, var = "FERTILIZER")

The above can be useful for exploring the full posteriors of all the parameters of performing specific calculations using the posteriors, summarising the parameters takes a few more steps.

fert.brm3 |>
  tidy_draws() |>
  gather_variables() |>
  median_hdci()
# A tibble: 15 × 7
   .variable                     .value   .lower  .upper .width .point .interval
   <chr>                          <dbl>    <dbl>   <dbl>  <dbl> <chr>  <chr>    
 1 accept_stat__                 0.961    0.712    1       0.95 median hdci     
 2 b_Intercept                 163.     150.     179.      0.95 median hdci     
 3 b_scaleFERTILIZERscaleEQFA…   0.806    0.599    1.00    0.95 median hdci     
 4 divergent__                   0        0        0       0.95 median hdci     
 5 energy__                     54.6     52.0     58.5     0.95 median hdci     
 6 Intercept                   163.     150.     179.      0.95 median hdci     
 7 lp__                        -53.0    -56.1    -51.6     0.95 median hdci     
 8 lprior                      -12.0    -12.1    -11.9     0.95 median hdci     
 9 n_leapfrog__                  7        3        7       0.95 median hdci     
10 prior_b                      -0.0607  -3.91     3.96    0.95 median hdci     
11 prior_Intercept             161.     -10.8    333.      0.95 median hdci     
12 prior_sigma                  68.8      0.0198 282.      0.95 median hdci     
13 sigma                        21.3     12.4     36.2     0.95 median hdci     
14 stepsize__                    0.588    0.521    0.594   0.95 median hdci     
15 treedepth__                   3        2        3       0.95 median hdci     

The gather_draws on the other hand, conveniently combines tidy_draws and gather_variables together in a single command. Furthermore, it returns all of the variables. The spread_draws function allows users to target specific variables (either by naming them in full or via regexp).

fert.brm3 |> get_variables()
 [1] "b_Intercept"                   "b_scaleFERTILIZERscaleEQFALSE"
 [3] "sigma"                         "Intercept"                    
 [5] "prior_Intercept"               "prior_b"                      
 [7] "prior_sigma"                   "lprior"                       
 [9] "lp__"                          "accept_stat__"                
[11] "stepsize__"                    "treedepth__"                  
[13] "n_leapfrog__"                  "divergent__"                  
[15] "energy__"                     
fert.brm3 |>
  gather_draws(b_Intercept, b_scaleFERTILIZERscaleEQFALSE, sigma) |>
  median_hdci()
# A tibble: 3 × 7
  .variable                      .value  .lower .upper .width .point .interval
  <chr>                           <dbl>   <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 b_Intercept                   163.    150.    179.     0.95 median hdci     
2 b_scaleFERTILIZERscaleEQFALSE   0.806   0.599   1.00   0.95 median hdci     
3 sigma                          21.3    12.4    36.2    0.95 median hdci     
## OR via regex
fert.brm3 |>
  gather_draws(`b_.*|sigma`, regex = TRUE) |>
  median_hdci()
# A tibble: 3 × 7
  .variable                      .value  .lower .upper .width .point .interval
  <chr>                           <dbl>   <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 b_Intercept                   163.    150.    179.     0.95 median hdci     
2 b_scaleFERTILIZERscaleEQFALSE   0.806   0.599   1.00   0.95 median hdci     
3 sigma                          21.3    12.4    36.2    0.95 median hdci     
fert.brm3 |>
  gather_draws(b_Intercept, b_scaleFERTILIZERscaleEQFALSE, sigma) |>
  ggplot() +
  stat_halfeye(aes(x = .value, y = .variable)) +
  facet_wrap(~.variable, scales = "free")

fert.brm3 |>
  gather_draws(b_Intercept, b_scaleFERTILIZERscaleEQFALSE, sigma) |>
  ggplot() +
  stat_halfeye(aes(x = .value, y = .variable))

Conclusions:

  • the estimated intercept (expected yield when fertilizer concentration is 0) is 163.41. This is the median of the posterior distribution for this parameter.
  • the estimated slope (rate at which yield changes per 1 unit change in fertilizer concentration), is 0.81 (mean) or (median). The 95% credibility intervals indicate that we are 95% confident that the slope is between 0.6 and 1 - e.g. there is a significant positive trend.
  • sigma is estimated to be 21.34

Select just the parameters of interest.

fert.brm3 |> spread_draws(b_Intercept, b_scaleFERTILIZERscaleEQFALSE, sigma)
# A tibble: 2,400 × 6
   .chain .iteration .draw b_Intercept b_scaleFERTILIZERscaleEQFALSE sigma
    <int>      <int> <int>       <dbl>                         <dbl> <dbl>
 1      1          1     1        166.                         0.890  16.2
 2      1          2     2        167.                         0.954  23.2
 3      1          3     3        161.                         0.796  14.8
 4      1          4     4        162.                         0.826  16.1
 5      1          5     5        163.                         0.777  18.6
 6      1          6     6        161.                         0.778  14.7
 7      1          7     7        167.                         0.833  20.6
 8      1          8     8        163.                         0.871  26.8
 9      1          9     9        160.                         0.826  23.2
10      1         10    10        140.                         0.584  26.4
# ℹ 2,390 more rows
# OR via regex
fert.brm3 |> spread_draws(`b_.*|sigma`, regex = TRUE)
# A tibble: 2,400 × 6
   .chain .iteration .draw b_Intercept b_scaleFERTILIZERscaleEQFALSE sigma
    <int>      <int> <int>       <dbl>                         <dbl> <dbl>
 1      1          1     1        166.                         0.890  16.2
 2      1          2     2        167.                         0.954  23.2
 3      1          3     3        161.                         0.796  14.8
 4      1          4     4        162.                         0.826  16.1
 5      1          5     5        163.                         0.777  18.6
 6      1          6     6        161.                         0.778  14.7
 7      1          7     7        167.                         0.833  20.6
 8      1          8     8        163.                         0.871  26.8
 9      1          9     9        160.                         0.826  23.2
10      1         10    10        140.                         0.584  26.4
# ℹ 2,390 more rows
fert.brm3 |> mcmc_plot(type = "intervals")

Note, this is now deprecated

fert.brm3 |>
  posterior_samples() |>
  as_tibble()
Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.
# A tibble: 2,400 × 9
   b_Intercept b_scaleFERTILIZERscaleE…¹ sigma Intercept prior_Intercept prior_b
         <dbl>                     <dbl> <dbl>     <dbl>           <dbl>   <dbl>
 1        166.                     0.890  16.2      166.           213.  -2.09  
 2        167.                     0.954  23.2      167.           135.   1.38  
 3        161.                     0.796  14.8      161.           293.   1.67  
 4        162.                     0.826  16.1      162.           241.  -4.03  
 5        163.                     0.777  18.6      163.           152.   0.221 
 6        161.                     0.778  14.7      161.           133.   0.0342
 7        167.                     0.833  20.6      167.            76.9 -2.71  
 8        163.                     0.871  26.8      163.           144.   0.805 
 9        160.                     0.826  23.2      160.           107.  -1.15  
10        140.                     0.584  26.4      140.           220.   1.70  
# ℹ 2,390 more rows
# ℹ abbreviated name: ¹​b_scaleFERTILIZERscaleEQFALSE
# ℹ 3 more variables: prior_sigma <dbl>, lprior <dbl>, lp__ <dbl>
fert.brm3 |>
  bayes_R2(summary = FALSE) |>
  median_hdci()
          y      ymin      ymax .width .point .interval
1 0.9185581 0.8246274 0.9273001   0.95 median      hdci

It is possible to obtain an empirical p-value (for those that cannot live without p-values). This is essentially compares (for any posterior) that the mean is zero compared to a multivariate distribution with elliptical contours.

mcmcpvalue <- function(samp) {
  ## elementary version that creates an empirical p-value for the
  ## hypothesis that the columns of samp have mean zero versus a
  ## general multivariate distribution with elliptical contours.

  ## differences from the mean standardized by the observed
  ## variance-covariance factor

  ## Note, I put in the bit for single terms
  if (length(dim(samp)) == 0) {
    std <- backsolve(chol(var(samp)), cbind(0, t(samp)) - mean(samp), transpose = TRUE)
    sqdist <- colSums(std * std)
    sum(sqdist[-1] > sqdist[1]) / length(samp)
  } else {
    std <- backsolve(chol(var(samp)), cbind(0, t(samp)) - colMeans(samp), transpose = TRUE)
    sqdist <- colSums(std * std)
    sum(sqdist[-1] > sqdist[1]) / nrow(samp)
  }
}

mcmcpvalue(as.matrix(fert.brm3)[, c("b_scaleFERTILIZERscaleEQFALSE")])
[1] 0
fert.brm3 |> modelsummary(
  statistic = c("conf.low", "conf.high"),
  shape = term ~ statistic
)
Warning: 
`modelsummary` uses the `performance` package to extract goodness-of-fit
statistics from models of this class. You can specify the statistics you wish
to compute by supplying a `metrics` argument to `modelsummary`, which will then
push it forward to `performance`. Acceptable values are: "all", "common",
"none", or a character vector of metrics names. For example: `modelsummary(mod,
metrics = c("RMSE", "R2")` Note that some metrics are computationally
expensive. See `?performance::performance` for details.
 This warning appears once per session.
(1)
Est. 2.5 % 97.5 %
b_Intercept 163.414 148.725 178.636
b_scaleFERTILIZERscaleEQFALSE 0.806 0.604 1.013
sigma 21.343 13.411 39.225
Num.Obs. 10
R2 0.919
R2 Adj. 0.890
ELPD -45.8
ELPD s.e. 1.4
LOOIC 91.6
LOOIC s.e. 2.9
WAIC 91.3
RMSE 17.18
fert.brm3 |> modelplot()

INLA captures a summary of the parameter posteriors within the fitted object.

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

fert.inla1 |> summary()
Time used:
    Pre = 0.247, Running = 0.263, Post = 0.0297, Total = 0.54 
Fixed effects:
              mean     sd 0.025quant 0.5quant 0.975quant   mode kld
(Intercept) 52.017 13.486     25.049   52.007     79.046 52.009   0
FERTILIZER   0.811  0.087      0.637    0.811      0.985  0.811   0

Model hyperparameters:
                                         mean    sd 0.025quant 0.5quant
Precision for the Gaussian observations 0.003 0.001      0.001    0.003
                                        0.975quant  mode
Precision for the Gaussian observations      0.007 0.002

Deviance Information Criterion (DIC) ...............: 91.27
Deviance Information Criterion (DIC, saturated) ....: 13.40
Effective number of parameters .....................: 3.09

Watanabe-Akaike information criterion (WAIC) ...: 91.15
Effective number of parameters .................: 2.51

Marginal log-Likelihood:  -55.92 
CPO, PIT is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

Conclusions:

  • in the initial block of information, we are reminded of the function call
  • a breakdown of the various components of the time taken for the analysis to run
  • the ‘Fixed effects’ indicate that the:
    • estimated Intercept is 52.017 (mean) or 52.007 (median)
    • estimated slope is 0.811 (mean) or 0.811 (median)
    • the 95% Credibility intervals indicate that we are 95% confident that the slope is between 0.637 and 0.985
  • the `Model hyperparameters’ indicate that the:
    • the precision of the Gaussian distribution is 0.003 (mean) and 0.003 (median) . Precision is \(\tau = 1/\sigma^2\) and thus, \(\sigma\) is estimated to be 18.2574186
    • WAIC (Watanabe-Akaike information criterion is estimated to be 91.148
brinla::bri.fixed.plot(fert.inla1)

brinla::bri.hyperpar.plot(fert.inla1)

To focus on and filter to just the fixed effects, we could now:

fert.draws |>
  filter(grepl("Intercept|FERTILIZER", Parameter)) |>
  group_by(Parameter) |>
  median_hdci(value)
# A tibble: 2 × 7
  Parameter      value .lower .upper .width .point .interval
  <chr>          <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 (Intercept):1 52.3   24.0   79.3     0.95 median hdci     
2 FERTILIZER:1   0.806  0.643  0.989   0.95 median hdci     
fert.draws |>
  filter(grepl("Intercept|FERTILIZER", Parameter)) |>
  ggplot(aes(x = value, y = Parameter)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  stat_pointinterval(point_interval = median_hdci)

## or separate into panels
fert.draws |>
  filter(grepl("Intercept|FERTILIZER", Parameter)) |>
  ggplot(aes(x = value)) +
  geom_vline(xintercept = 0) +
  stat_pointinterval(point_interval = median_hdci) +
  facet_grid(~Parameter, scales = "free_x")

fert.draws |>
  filter(grepl("Intercept|FERTILIZER", Parameter)) |>
  group_by(Parameter) |>
  nest() |>
  mutate(
    hpd = map(data, ~ median_hdci(.x$value)),
    Density = map(data, function(.x) {
      d <- density(.x$value)
      data.frame(x = d$x, y = d$y)
    }),
    Quant = map2(Density, hpd, ~ factor(findInterval(.x$x, .y[, c("ymin", "ymax")])))
  ) |>
  unnest(c(Density, Quant)) |>
  ggplot(aes(x = x, y = y, fill = Quant)) +
  geom_vline(xintercept = 0) +
  geom_ribbon(aes(ymin = 0, ymax = y)) +
  facet_wrap(~Parameter, scales = "free")

fert.draws |>
  filter(grepl("Intercept|FERTILIZER", Parameter)) |>
  ggplot(aes(x = value, y = Parameter, fill = stat(quantile))) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  ggridges::stat_density_ridges(
    geom = "density_ridges_gradient",
    calc_ecdf = TRUE,
    quantile_fun = function(x, probs) quantile(x, probs),
    quantiles = c(0.025, 0.975)
  ) +
  scale_fill_viridis_d()
Warning: `stat(quantile)` was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(quantile)` instead.
Picking joint bandwidth of 1.43

fert.draws |>
  filter(grepl("Intercept|FERTILIZER", Parameter)) |>
  ggplot(aes(x = value, y = Parameter, fill = stat(quantile))) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  ggridges::stat_density_ridges(
    geom = "density_ridges_gradient",
    calc_ecdf = TRUE,
    quantile_fun = function(x, probs) quantile(x, probs),
    quantiles = c(0.025, 0.975)
  ) +
  scale_fill_viridis_d() +
  facet_grid(~Parameter, scales = "free_x")
Picking joint bandwidth of 2.85
Picking joint bandwidth of 0.0191

cor(
  fert.inla1$summary.fitted.values[, "mean"],
  fert$YIELD
)^2
[1] 0.9216005

10 Predictions

There are a large number of candidate routines for performing prediction. We will go through many of these. It is worth noting that prediction is technically the act of estimating what we expect to get if we were to collect a single new observation from a particular population (e.g. a specific level of fertilizer concentration). Often this is not what we want. Often we want the fitted values - estimates of what we expect to get if we were to collect multiple new observations and average them.

So while fitted values represent the expected underlying processes occurring in the system, predicted values represent our expectations from sampling from such processes.

Package Function Description Summarise with
rstantools posterior_predict Draw from the posterior of a prediction (includes sigma) - predicts single observations tidyMCMC()
rstantools posterior_linpred Draw from the posterior of the fitted values (on the link scale) - predicts average observations tidyMCMC()
rstantools posterior_epred Draw from the posterior of the fitted values (on the response scale) - predicts average observations tidyMCMC()
tidybayes predicted_draws Extract the posterior of prediction values median_hdci()
tidybayes epred_draws Extract the posterior of expected values median_hdci()
tidybayes fitted_draws Extract the posterior of fitted values median_hdci()
tidybayes add_predicted_draws Adds draws from the posterior of predictions to a data frame (of prediction data) median_hdci()
tidybayes add_fitted_draws Adds draws from the posterior of fitted values to a data frame (of prediction data) median_hdci()
emmeans emmeans Estimated marginal means from which posteriors can be drawn (via tidy_draws median_hdci()

For simple models prediction is essentially taking the model formula complete with parameter (coefficient) estimates and solving for new values of the predictor. To explore this, we will use the fitted model to predict Yield for a Fertilizer concentration of 110.

We will therefore start by establishing this prediction domain as a data frame to use across all of the prediction routines.

## establish a data set that defines the new data to predict against
newdata <- data.frame(FERTILIZER = 110)
fert.rstanarm3 |> emmeans(~FERTILIZER, at = newdata)
 FERTILIZER emmean lower.HPD upper.HPD
        110    141       127       157

Point estimate displayed: median 
HPD interval probability: 0.95 
fert.rstanarm3 |>
  tidybayes::predicted_draws(newdata) |>
  median_hdci()
# A tibble: 1 × 8
  FERTILIZER  .row .prediction .lower .upper .width .point .interval
       <dbl> <int>       <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1        110     1        141.   86.4   189.   0.95 median hdci     
## or
newdata |>
  tidybayes::add_predicted_draws(fert.rstanarm3) |>
  median_hdci()
# A tibble: 1 × 8
  FERTILIZER  .row .prediction .lower .upper .width .point .interval
       <dbl> <int>       <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1        110     1        141.   93.3   192.   0.95 median hdci     
## or
fert.rstanarm3 |>
  posterior_predict(newdata = newdata) |>
  tidyMCMC(estimate.method = "median", conf.int = TRUE, conf.method = "HPDinterval")
# A tibble: 1 × 5
  term  estimate std.error conf.low conf.high
  <chr>    <dbl>     <dbl>    <dbl>     <dbl>
1 1         142.      24.5     94.7      191.
fert.rstanarm3 |>
  epred_draws(newdata) |>
  median_hdci()
# A tibble: 1 × 8
  FERTILIZER  .row .epred .lower .upper .width .point .interval
       <dbl> <int>  <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1        110     1   141.   127.   157.   0.95 median hdci     
## or
newdata |>
  add_epred_draws(fert.rstanarm3) |>
  median_hdci()
# A tibble: 1 × 8
  FERTILIZER  .row .epred .lower .upper .width .point .interval
       <dbl> <int>  <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1        110     1   141.   127.   157.   0.95 median hdci     
## or
fert.rstanarm3 |>
  posterior_epred(newdata = newdata) |>
  tidyMCMC(estimate.method = "median", conf.int = TRUE, conf.method = "HPDinterval")
# A tibble: 1 × 5
  term  estimate std.error conf.low conf.high
  <chr>    <dbl>     <dbl>    <dbl>     <dbl>
1 1         141.      7.75     127.      157.
fert.rstanarm3 |>
  linpred_draws(newdata) |>
  median_hdci()
# A tibble: 1 × 8
  FERTILIZER  .row .linpred .lower .upper .width .point .interval
       <dbl> <int>    <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1        110     1     141.   127.   157.   0.95 median hdci     
## or
newdata |>
  add_linpred_draws(fert.rstanarm3) |>
  median_hdci()
# A tibble: 1 × 8
  FERTILIZER  .row .linpred .lower .upper .width .point .interval
       <dbl> <int>    <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1        110     1     141.   127.   157.   0.95 median hdci     
## or
fert.rstanarm3 |>
  posterior_linpred(newdata = newdata) |>
  tidyMCMC(estimate.method = "median", conf.int = TRUE, conf.method = "HPDinterval")
# A tibble: 1 × 5
  term  estimate std.error conf.low conf.high
  <chr>    <dbl>     <dbl>    <dbl>     <dbl>
1 1         141.      7.75     127.      157.
## Establish an appropriate model matrix
Xmat <- model.matrix(~FERTILIZER, data = newdata)
### get the posterior draws for the linear predictor
coefs <- fert.rstanarm3$stanfit |>
  as_draws_df() |>
  dplyr::select(c("(Intercept)", "FERTILIZER")) |>
  as.matrix()
Warning: Dropping 'draws_df' class as required metadata was removed.
fit <- coefs %*% t(Xmat)
fit |> median_hdci()
         y     ymin     ymax .width .point .interval
1 141.1377 127.4528 157.2793   0.95 median      hdci
## or
coefs <- fert.rstanarm3 |>
  tidy_draws() |>
  dplyr::select(c("(Intercept)", "FERTILIZER")) |>
  as.matrix()
fit <- coefs %*% t(Xmat)
fit |> median_hdci()
         y     ymin     ymax .width .point .interval
1 141.1377 127.4528 157.2793   0.95 median      hdci
## or
coefs <- fert.rstanarm3$stanfit |>
  as_draws_matrix() |>
  subset_draws(c("(Intercept)", "FERTILIZER"))
fit <- coefs %*% t(Xmat)
fit |> median_hdci()
         y     ymin     ymax .width .point .interval
1 141.1377 127.4528 157.2793   0.95 median      hdci

This is the option we will focus on for most of the worksheets

fert.brm3 |>
  emmeans(~FERTILIZER, at = newdata)
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
 FERTILIZER emmean lower.HPD upper.HPD
        110    141       126       158

Point estimate displayed: median 
HPD interval probability: 0.95 
## OR
fert.brm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  tidy_draws() |>
  median_hdci()
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
# A tibble: 1 × 6
  `FERTILIZER 110` .lower .upper .width .point .interval
             <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1             141.   126.   158.   0.95 median hdci     
fert.mcmc <- fert.brm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  tidy_draws()
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

What is the probability that if we apply a fertiliser concentration of 110 units, the yield will exceed 120 units?

fert.brm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  tidy_draws() |>
  ggplot(aes(x = `FERTILIZER 110`)) +
  geom_vline(xintercept = 120) +
  geom_density(fill = "orange", alpha = 0.3)
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

fert.brm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  tidy_draws() |>
  summarise_draws(
    median,
    HDInterval::hdi,
    Pg = ~ mean(.x > 120)
  )
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
# A tibble: 1 × 5
  variable       median lower upper    Pg
  <chr>           <dbl> <dbl> <dbl> <dbl>
1 FERTILIZER 110   141.  126.  158. 0.992
fert.brm3 |>
  tidybayes::predicted_draws(newdata) |>
  median_hdci()
# A tibble: 1 × 8
  FERTILIZER  .row .prediction .lower .upper .width .point .interval
       <dbl> <int>       <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1        110     1        141.   95.0   195.   0.95 median hdci     
## or
newdata |>
  tidybayes::add_predicted_draws(fert.brm3) |>
  median_hdci()
# A tibble: 1 × 8
  FERTILIZER  .row .prediction .lower .upper .width .point .interval
       <dbl> <int>       <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1        110     1        141.   95.2   194.   0.95 median hdci     
## or
fert.brm3 |>
  posterior_predict(newdata = newdata) |>
  as.mcmc() |>
  tidyMCMC(estimate.method = "median", conf.int = TRUE, conf.method = "HPDinterval")
# A tibble: 1 × 5
  term  estimate std.error conf.low conf.high
  <chr>    <dbl>     <dbl>    <dbl>     <dbl>
1 var1      141.      24.4     91.1      188.
fert.brm3 |>
  epred_draws(newdata) |>
  median_hdci()
# A tibble: 1 × 8
  FERTILIZER  .row .epred .lower .upper .width .point .interval
       <dbl> <int>  <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1        110     1   141.   126.   158.   0.95 median hdci     
## or
newdata |>
  add_epred_draws(fert.brm3) |>
  median_hdci()
# A tibble: 1 × 8
  FERTILIZER  .row .epred .lower .upper .width .point .interval
       <dbl> <int>  <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1        110     1   141.   126.   158.   0.95 median hdci     
## or
fert.brm3 |>
  posterior_epred(newdata = newdata) |>
  as.mcmc() |>
  tidyMCMC(estimate.method = "median", conf.int = TRUE, conf.method = "HPDinterval")
# A tibble: 1 × 5
  term  estimate std.error conf.low conf.high
  <chr>    <dbl>     <dbl>    <dbl>     <dbl>
1 var1      141.      7.94     126.      158.
fert.brm3 |>
  linpred_draws(newdata) |>
  median_hdci()
# A tibble: 1 × 8
  FERTILIZER  .row .linpred .lower .upper .width .point .interval
       <dbl> <int>    <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1        110     1     141.   126.   158.   0.95 median hdci     
## or
newdata |>
  add_linpred_draws(fert.brm3) |>
  median_hdci()
# A tibble: 1 × 8
  FERTILIZER  .row .linpred .lower .upper .width .point .interval
       <dbl> <int>    <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1        110     1     141.   126.   158.   0.95 median hdci     
## or
fert.brm3 |>
  posterior_linpred(newdata = newdata) |>
  as.mcmc() |>
  tidyMCMC(estimate.method = "median", conf.int = TRUE, conf.method = "HPDinterval")
# A tibble: 1 × 5
  term  estimate std.error conf.low conf.high
  <chr>    <dbl>     <dbl>    <dbl>     <dbl>
1 var1      141.      7.94     126.      158.
## Establish an appropriate model matrix
Xmat <- model.matrix(~FERTILIZER, data = newdata)
### get the posterior draws for the linear predictor
coefs <- fert.brm3 |>
  as_draws_df() |>
  dplyr::select(c("b_Intercept", "b_scaleFERTILIZERscaleEQFALSE")) |>
  as.matrix()
Warning: Dropping 'draws_df' class as required metadata was removed.
fit <- coefs %*% t(Xmat)
fit |>
  median_hdci() |>
  bind_cols(newdata)
         y     ymin    ymax .width .point .interval FERTILIZER
1 251.9218 224.2932 277.191   0.95 median      hdci        110
## or
coefs <- fert.brm3 |>
  tidy_draws() |>
  dplyr::select(c("b_Intercept", "b_scaleFERTILIZERscaleEQFALSE")) |>
  as.matrix()
fit <- coefs %*% t(Xmat)
fit |>
  median_hdci() |>
  bind_cols(newdata)
         y     ymin    ymax .width .point .interval FERTILIZER
1 251.9218 224.2932 277.191   0.95 median      hdci        110
## or
coefs <- fert.brm3 |>
  as_draws_matrix() |>
  subset_draws(c("b_Intercept", "b_scaleFERTILIZERscaleEQFALSE"))
fit <- coefs %*% t(Xmat)
fit |>
  median_hdci() |>
  bind_cols(newdata)
         y     ymin    ymax .width .point .interval FERTILIZER
1 251.9218 224.2932 277.191   0.95 median      hdci        110
## or
## Establish an appropriate model matrix
Xmat <- model.matrix(~FERTILIZER, data = newdata)
### get the posterior draws for the linear predictor
coefs <- fert.brm3 |>
  posterior_samples(pars = c("b_Intercept", "b_scaleFERTILIZERscaleEQFALSE")) |>
  as.matrix()
Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.
fit <- coefs %*% t(Xmat)
fit |>
  median_hdci() |>
  bind_cols(newdata)
         y     ymin    ymax .width .point .interval FERTILIZER
1 251.9218 224.2932 277.191   0.95 median      hdci        110
## Expected values
Xmat <- model.matrix(~FERTILIZER, newdata)
nms <- colnames(fert.inla1$model.matrix)
sel <- sapply(nms, function(x) 0, simplify = FALSE)
n <- 1000
draws <- inla.posterior.sample(n = n, result = fert.inla1, selection = sel, seed = 1)
Warning in inla.posterior.sample(n, rfake, intern = intern, use.improved.mean =
use.improved.mean, : Since 'seed!=0', parallel model is disabled and serial
model is selected, num.threads='1:1'
coefs <- t(sapply(draws, function(x) x$latent))
Fit <- coefs %*% t(Xmat) |>
  as.data.frame() |>
  mutate(Rep = 1:n()) |>
  pivot_longer(cols = -Rep) |>
  group_by(name) |>
  median_hdci(value) |>
  ungroup() |>
  mutate(name = as.integer(as.character(name))) |>
  arrange(name)
newdata.inla <- newdata |>
  bind_cols(Fit)
newdata.inla
  FERTILIZER name    value   .lower   .upper .width .point .interval
1        110    1 141.3302 129.1472 155.7588   0.95 median      hdci
## Predictions
Fit <- coefs %*% t(Xmat)
draws <- inla.posterior.sample(n = n, result = fert.inla1, seed = 1)
Warning in inla.posterior.sample(n, rfake, intern = intern, use.improved.mean =
use.improved.mean, : Since 'seed!=0', parallel model is disabled and serial
model is selected, num.threads='1:1'
sigma <- sqrt(1 / (sapply(draws, function(x) x$hyperpar)))
sigma <- sapply(sigma, function(x) rnorm(1, 0, sigma))
Fit <- sweep(Fit, MARGIN = 1, sigma, FUN = "+") |>
  as.data.frame() |>
  mutate(Rep = 1:n()) |>
  pivot_longer(cols = -Rep) |>
  group_by(name) |>
  median_hdci(value) |>
  bind_cols(newdata) |>
  ungroup() |>
  dplyr::select(FERTILIZER, everything(), -name)
Fit
# A tibble: 1 × 7
  FERTILIZER value .lower .upper .width .point .interval
       <dbl> <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1        110  142.   61.0   226.   0.95 median hdci     
## or
fun <- function(coefs = NA) {
  ## theta[1] is the precision
  return(Intercept + FERTILIZER * coefs[, "FERTILIZER"] +
    rnorm(nrow(coefs), sd = sqrt(1 / theta[1])))
}
Fit <- inla.posterior.sample.eval(fun, draws, coefs = newdata) |>
  as.data.frame() |>
  bind_cols(newdata) |>
  pivot_longer(cols = -FERTILIZER) |>
  group_by(FERTILIZER) |>
  median_hdci(value)

Fit
# A tibble: 1 × 7
  FERTILIZER value .lower .upper .width .point .interval
       <dbl> <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1        110  141.   101.   185.   0.95 median hdci     
fert.pred <- fert |>
  bind_rows(newdata)
i.newdata <- (nrow(fert) + 1):nrow(fert.pred)
fert.inla2 <- inla(YIELD ~ FERTILIZER,
  data = fert.pred,
  family = "gaussian",
  control.fixed = list(
    mean.intercept = 80,
    prec.intercept = 0.00001,
    mean = 0,
    prec = 0.0384
  ),
  control.family = list(hyper = list(prec = list(
    prior = "loggamma",
    param = c(0.5, 0.31)
  ))),
  control.compute = list(config = TRUE, dic = TRUE, waic = TRUE, cpo = TRUE)
)

fert.inla2$summary.fitted.values[i.newdata, ]
                        mean       sd 0.025quant 0.5quant 0.975quant    mode
fitted.Predictor.11 141.2114 6.598545   128.0136  141.209   154.4243 141.205
## or on the link scale...
fert.inla2$summary.linear.predictor[i.newdata, ]
                 mean       sd 0.025quant 0.5quant 0.975quant     mode
Predictor.11 141.2114 6.598567   128.0136  141.209   154.4243 141.2093
                      kld
Predictor.11 5.304516e-08
Xmat <- model.matrix(~FERTILIZER, data = newdata)
lincomb <- inla.make.lincombs(as.data.frame(Xmat))

fert.inla3 <- inla(YIELD ~ FERTILIZER,
  data = fert,
  family = "gaussian",
  lincomb = lincomb,
  control.fixed = list(
    mean.intercept = 80,
    prec.intercept = 0.00001,
    mean = 0,
    prec = 0.0384
  ),
  control.family = list(hyper = list(prec = list(
    prior = "loggamma",
    param = c(0.5, 0.31)
  ))),
  control.compute = list(config = TRUE, dic = TRUE, waic = TRUE, cpo = TRUE)
)

fert.inla3$summary.lincomb.derived
             ID     mean       sd 0.025quant 0.5quant 0.975quant     mode
lc0000000001  1 141.2114 6.598567   128.0136  141.209   154.4243 141.2093
                      kld
lc0000000001 5.304516e-08

11 Hypothesis testing

Since we have the entire posterior, we are able to make probability statements. We simply count up the number of MCMC sample draws that satisfy a condition (e.g represent a slope greater than 0) and then divide by the total number of MCMC samples.

For this exercise, we will explore the following:

  • what this the probability that Yield increases with increasing Fertilizer concentration - this is just asking what proportion of the drawn slopes are greater than 0.
  • how much of a change in yield do we expect if fertilizer concentration is increased from 100 to 200 units and what is the probability that the change is more than 50%
fert.rstanarm3 |> hypothesis("FERTILIZER>0")
Hypothesis Tests for class :
        Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
1 (FERTILIZER) > 0     0.81       0.1     0.65     0.97        Inf         1
  Star
1    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

Conclusions:

  • the probability of a positive effect of Fertilizer on grass Yield is 1. That is, there is, we are 100% confident that there is an effect.
  • the evidence ratio is normally the ratio of the number of cases that satisfy the hypothesis to the number of cases that do not. Since the number of cases that do not satisfy the hypothesis is 0, the evidence ratio is Inf (since division by 0)

Alternatively…

fert.rstanarm3 |>
  tidy_draws() |>
  summarise(P = mean(FERTILIZER > 0))
# A tibble: 1 × 1
      P
  <dbl>
1     1

Conclusions:

  • the probability of a positive effect of Fertilizer on grass Yield is 1. That is, there is, we are 100% confident that there is an effect.

The procedure highlighted above for calculating excedence probabilities evaluates the degree of evidence for a effect in a particular direction. However, there are other instances where there is a desire to evaluate the evidence that something has change (either increased or decreased). Such purposes are similar to the Frequentist pursuit of testing a null hypothesis (e.g. effect = 0).

The Region of Practical Equivalence (ROPE) evaluates evidence that an effect is “practically equivalent” to a value (e.g. 0) by calculating the proportion of effects that are within a nominated range. Kruschke (2018) argued that for standardized parameters, the range of -0.1 to 0.1 would envelop a negligible effect based on Cohen (1988). Kruschke (2018) also suggested that this range could be extended to non-standardized parameters by multiplying by the standard deviation of the response. Accordingly, calculating the proportion of posterior density within this ROPE could act as a form of “null-hypothesis” testing in a Bayesian framework.

0.1 * sd(fert$YIELD) / sd(fert$FERTILIZER)
[1] 0.08452018
## Cannot pipe to rope if want to pipe rope to plot()
bayestestR::rope(fert.rstanarm3, parameters = "FERTILIZER", range = c(-0.08, 0.08), ci_method = "HDI")
# Proportion of samples inside the ROPE [-0.08, 0.08]:

Parameter  | inside ROPE
------------------------
FERTILIZER |      0.00 %
bayestestR::rope(fert.rstanarm3, parameters = "FERTILIZER", range = c(-0.08, 0.08), ci_method = "HDI") |> plot()

bayestestR::equivalence_test(fert.rstanarm3, parameters = "FERTILIZER", range = c(-0.08, 0.08), ci_method = "HDI")
# Test for Practical Equivalence

  ROPE: [-0.08 0.08]

Parameter  |       H0 | inside ROPE |      95% HDI
--------------------------------------------------
FERTILIZER | Rejected |      0.00 % | [0.61, 1.01]
bayestestR::equivalence_test(fert.rstanarm3, parameters = "FERTILIZER", range = c(-0.08, 0.08), ci_method = "HDI") |> plot()
Picking joint bandwidth of 0.0162

newdata <- list(FERTILIZER = c(200, 100))
fert.rstanarm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  pairs()
 contrast                      estimate lower.HPD upper.HPD
 FERTILIZER200 - FERTILIZER100     80.9      61.6       102

Point estimate displayed: median 
HPD interval probability: 0.95 

Conclusions:

  • the change in Yield associated with an increase in Fertilizer concentration from 100 to 200 is 80.89.

We can alternatively express this as a percentage change…

newdata <- list(FERTILIZER = c(200, 100))
fert.rstanarm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  regrid(transform = "log") |>
  pairs() |>
  regrid()
 contrast                    ratio lower.HPD upper.HPD
 FERTILIZER200/FERTILIZER100  1.61      1.41      1.81

Point estimate displayed: median 
HPD interval probability: 0.95 

If we want to derive other properties, such as the percentage change, then we use tidy_draws() and then simple tidyverse spreadsheet operation.

fert.rstanarm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  regrid(transform = "log") |>
  pairs() |>
  regrid() |>
  tidy_draws() |>
  summarise_draws(median,
    HDInterval::hdi,
    P = ~ mean(.x > 1.5)
  )
# A tibble: 1 × 5
  variable                             median lower upper     P
  <chr>                                 <dbl> <dbl> <dbl> <dbl>
1 contrast FERTILIZER200/FERTILIZER100   1.61  1.41  1.81 0.884
fert.mcmc <- fert.rstanarm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  tidy_draws() |>
  rename_with(~ str_replace(., "FERTILIZER ", "p")) |>
  mutate(
    Eff = p200 - p100,
    PEff = 100 * Eff / p100
  )
fert.mcmc |> head()
# A tibble: 6 × 7
  .chain .iteration .draw  p200  p100   Eff  PEff
   <int>      <int> <int> <dbl> <dbl> <dbl> <dbl>
1      1          1     1  211.  125.  86.5  69.3
2      1          2     2  213.  131.  81.8  62.4
3      1          3     3  217.  124.  92.9  75.0
4      1          4     4  207.  130.  77.2  59.5
5      1          5     5  208.  125.  83.0  66.3
6      1          6     6  213.  121.  91.8  75.9

Now we can calculate the medians and HPD intervals of each column (and ignore the .chain, .iteration and .draw).

fert.mcmc |> tidyMCMC(
  estimate.method = "median",
  conf.int = TRUE, conf.method = "HPDinterval"
)
# A tibble: 7 × 5
  term       estimate std.error conf.low conf.high
  <chr>         <dbl>     <dbl>    <dbl>     <dbl>
1 .chain          2       0.817      1         3  
2 .iteration    400.    231.         1       761  
3 .draw        1200.    693.         1      2281  
4 p200          214.      9.62     196.      234. 
5 p100          133.      8.17     118.      149. 
6 Eff            80.8    10.1       61.6     102. 
7 PEff           61.2    10.00      41.4      80.8

Alternatively, we could use median_hdci

fert.mcmc |> median_hdci(PEff)
# A tibble: 1 × 6
   PEff .lower .upper .width .point .interval
  <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1  61.0   41.6   81.0   0.95 median hdci     
fert.mcmc |> median_hdci(PEff, Eff)
# A tibble: 1 × 9
   PEff PEff.lower PEff.upper   Eff Eff.lower Eff.upper .width .point .interval
  <dbl>      <dbl>      <dbl> <dbl>     <dbl>     <dbl>  <dbl> <chr>  <chr>    
1  61.0       41.6       81.0  80.9      61.9      102.   0.95 median hdci     

Conclusions:

  • the percentage change in Yield associated with an increase in Fertilizer concentration from 100 to 200 is 60.96% with a 95% credibility interval of 41.61 81

To get the probability that the effect is greater than a 50% increase.

fert.mcmc |> summarise(P = mean(PEff > 50))
# A tibble: 1 × 1
      P
  <dbl>
1 0.884

Conclusions:

  • the probability that Yield will increase by more than 50 following an increase in Fertilizer concentration from 100 to 200 is 0.88.

Finally, we could alternatively use hypothesis(). Note that when we do so, the estimate is the difference between the effect and the hypothesised effect (50%).

fert.mcmc |> hypothesis("PEff>50")
Hypothesis Tests for class :
       Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
1 (PEff)-(50) > 0    11.18        10    -4.34    27.86       7.63      0.88
  Star
1     
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
newdata <- list(FERTILIZER = c(200, 100))
fert.rstanarm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  pairs() |>
  tidy_draws() |>
  summarise(across(contains("contrast"),
    list(
      P = ~ mean(. > 50),
      HDCI = ~ median_hdci(.)
    ),
    .names = c("{.fn}")
  )) |>
  tidyr::unpack(HDCI)
# A tibble: 1 × 7
      P     y  ymin  ymax .width .point .interval
  <dbl> <dbl> <dbl> <dbl>  <dbl> <chr>  <chr>    
1 0.996  80.9  61.9  102.   0.95 median hdci     

We can also express this as a percentage change

newdata <- list(FERTILIZER = c(200, 100))
## Simple
fert.rstanarm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  regrid(transform = "log") |>
  pairs() |>
  regrid()
 contrast                    ratio lower.HPD upper.HPD
 FERTILIZER200/FERTILIZER100  1.61      1.41      1.81

Point estimate displayed: median 
HPD interval probability: 0.95 
## More advanced (both P and percent change)
fert.mcmc <- fert.rstanarm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  regrid(transform = "log") |>
  pairs() |>
  regrid() |>
  tidy_draws() |>
  mutate(across(contains("contrast"), ~ 100 * (. - 1)))

fert.mcmc |>
  summarise(across(contains("contrast"),
    list(
      P = ~ mean(. > 50),
      HDCI = ~ median_hdci(.)
    ),
    .names = c("{.fn}")
  )) |>
  tidyr::unpack(HDCI)
# A tibble: 1 × 7
      P     y  ymin  ymax .width .point .interval
  <dbl> <dbl> <dbl> <dbl>  <dbl> <chr>  <chr>    
1 0.884  61.0  41.6  81.0   0.95 median hdci     
## OR
fert.mcmc |>
  summarise_draws(median,
    HDInterval::hdi,
    P = ~ mean(.x > 50)
  )
# A tibble: 1 × 5
  variable                             median lower upper     P
  <chr>                                 <dbl> <dbl> <dbl> <dbl>
1 contrast FERTILIZER200/FERTILIZER100   61.0  41.4  80.8 0.884
fert.rstanarm3 |>
  linpred_draws(newdata = as.data.frame(newdata)) |>
  ungroup() |>
  group_by(.draw) |>
  summarise(
    Eff = -1 * diff(.linpred),
    PEff = 100 * Eff / .linpred[2]
  ) |>
  ungroup() |>
  mutate(P = mean(PEff > 50)) |>
  pivot_longer(cols = -.draw) |>
  group_by(name) |>
  median_hdci()
# A tibble: 3 × 7
  name   value .lower  .upper .width .point .interval
  <chr>  <dbl>  <dbl>   <dbl>  <dbl> <chr>  <chr>    
1 Eff   80.9   61.9   102.      0.95 median hdci     
2 P      0.884  0.884   0.884   0.95 median hdci     
3 PEff  61.0   41.6    81.0     0.95 median hdci     
## OR
fert.rstanarm3 |>
  epred_draws(newdata = as.data.frame(newdata)) |>
  ungroup() |>
  group_by(.draw) |>
  summarise(
    Eff = -1 * diff(.epred),
    PEff = 100 * Eff / .epred[2]
  ) |>
  ungroup() |>
  mutate(P = mean(PEff > 50)) |>
  pivot_longer(cols = -.draw) |>
  group_by(name) |>
  median_hdci()
# A tibble: 3 × 7
  name   value .lower  .upper .width .point .interval
  <chr>  <dbl>  <dbl>   <dbl>  <dbl> <chr>  <chr>    
1 Eff   80.9   61.9   102.      0.95 median hdci     
2 P      0.884  0.884   0.884   0.95 median hdci     
3 PEff  61.0   41.6    81.0     0.95 median hdci     
## OR for prediction of individual values
fert.rstanarm3 |>
  predicted_draws(newdata = as.data.frame(newdata)) |>
  ungroup() |>
  group_by(.draw) |>
  summarise(
    Eff = -1 * diff(.prediction),
    PEff = 100 * Eff / .prediction[2]
  ) |>
  ungroup() |>
  mutate(P = mean(PEff > 50)) |>
  pivot_longer(cols = -.draw) |>
  group_by(name) |>
  median_hdci()
# A tibble: 3 × 7
  name   value .lower  .upper .width .point .interval
  <chr>  <dbl>  <dbl>   <dbl>  <dbl> <chr>  <chr>    
1 Eff   81.3   11.7   150.      0.95 median hdci     
2 P      0.638  0.638   0.638   0.95 median hdci     
3 PEff  61.2    2.80  156.      0.95 median hdci     
fert.brm3 |> hypothesis("scaleFERTILIZERscaleEQFALSE > 0")
Hypothesis Tests for class b:
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (scaleFERTILIZERs... > 0     0.81       0.1     0.64     0.97        Inf
  Post.Prob Star
1         1    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
fert.brm3 |>
  hypothesis("scaleFERTILIZERscaleEQFALSE > 0") |>
  plot(ignore_prior = TRUE)

fert.brm3 |>
  hypothesis("scaleFERTILIZERscaleEQFALSE > 0") |>
  _[["samples"]] |>
  median_hdci()
# A tibble: 1 × 6
     H1 .lower .upper .width .point .interval
  <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 0.806  0.599   1.00   0.95 median hdci     

Conclusions:

  • the probability of a positive effect Fertilizer on grass yield is 1. That is, there is, we are 100% confident that there is an effect.
  • the evidence ratio is normally the ratio of the number of cases that satisfy the hypothesis to the number of cases that do not. Since the number of cases that do not satisfy the hypothesis is 0, the evidence ratio is Inf (since division by 0)

Alternatively…

fert.brm3 |> get_variables()
 [1] "b_Intercept"                   "b_scaleFERTILIZERscaleEQFALSE"
 [3] "sigma"                         "Intercept"                    
 [5] "prior_Intercept"               "prior_b"                      
 [7] "prior_sigma"                   "lprior"                       
 [9] "lp__"                          "accept_stat__"                
[11] "stepsize__"                    "treedepth__"                  
[13] "n_leapfrog__"                  "divergent__"                  
[15] "energy__"                     
fert.brm3 |>
  gather_draws(b_scaleFERTILIZERscaleEQFALSE) |>
  summarise(P = mean(.value > 0))
# A tibble: 1 × 2
  .variable                         P
  <chr>                         <dbl>
1 b_scaleFERTILIZERscaleEQFALSE     1
# OR
fert.brm3 |>
  tidy_draws() |>
  summarise(P = mean(b_scaleFERTILIZERscaleEQFALSE > 0))
# A tibble: 1 × 1
      P
  <dbl>
1     1
# OR
fert.brm3 |>
  tidy_draws() |>
  ## mutate(FERTILIZER = b_scaleFERTILIZER / fert_sigma) |>
  ## dplyr::select(matches("^FERT")) |>
  dplyr::select(matches("^b_scaleFERT")) |>
  summarise_draws(median,
    P = ~ mean(.x > 0)
  )
# A tibble: 1 × 3
  variable                      median     P
  <chr>                          <dbl> <dbl>
1 b_scaleFERTILIZERscaleEQFALSE  0.806     1

Conclusions:

  • the probability of a positive effect of Fertilizer on grass Yield is 1. That is, there is, we are 100% confident that there is an effect.
fert.brm3 |> hypothesis("scaleFERTILIZERscaleEQFALSE > 0.9")
Hypothesis Tests for class b:
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (scaleFERTILIZERs... > 0    -0.09       0.1    -0.26     0.07       0.18
  Post.Prob Star
1      0.15     
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
fert.brm3 |>
  hypothesis("scaleFERTILIZERscaleEQFALSE > 0.9") |>
  plot(ignore_prior = TRUE)

# This returns a list
fert.brm3 |>
  hypothesis("scaleFERTILIZERscaleEQFALSE > 0.9") |>
  plot(ignore_prior = TRUE) |>
  _[[1]] +
  geom_vline(xintercept = 0, linetype = "dashed")

# OR
fert.brm3 |>
  tidy_draws() |>
  ## mutate(FERTILIZER = b_scaleFERTILIZER / fert_sigma) |>
  ## dplyr::select(matches("^FERT")) |>
  dplyr::select(matches("^b_scaleFERT")) |>
  summarise_draws(median,
    P = ~ mean(.x > 0.7)
  )
# A tibble: 1 × 3
  variable                      median     P
  <chr>                          <dbl> <dbl>
1 b_scaleFERTILIZERscaleEQFALSE  0.806  0.87
fert.brm3 |>
  tidy_draws() |>
  ## mutate(FERTILIZER = b_scaleFERTILIZER / fert_sigma) |>
  mutate(FERTILIZER = b_scaleFERTILIZERscaleEQFALSE) |>
  ggplot(aes(x = FERTILIZER)) +
  geom_density(fill = "orange") +
  geom_vline(xintercept = 0.7, linetype = "dashed")

The procedure highlighted above for calculating excedence probabilities evaluates the degree of evidence for a effect in a particular direction. However, there are other instances where there is a desire to evaluate the evidence that something has change (either increased or decreased). Such purposes are similar to the Frequentist pursuit of testing a null hypothesis (e.g. effect = 0).

The Region of Practical Equivalence (ROPE) evaluates evidence that an effect is “practically equivalent” to a value (e.g. 0) by calculating the proportion of effects that are within a nominated range. Kruschke (2018) argued that for standardized parameters, the range of -0.1 to 0.1 would envelop a negligible effect based on Cohen (1988). Kruschke (2018) also suggested that this range could be extended to non-standardized parameters by multiplying by the standard deviation of the response. Accordingly, calculating the proportion of posterior density within this ROPE could act as a form of “null-hypothesis” testing in a Bayesian framework.

0.1 * sd(fert$YIELD) / sd(fert$FERTILIZER)
[1] 0.08452018
0.1 * sd(fert$YIELD)
[1] 6.397439
## Cannot pipe to rope if want to pipe rope to plot()
bayestestR::rope(fert.brm3,
  parameters = "FERTILIZER",
  range = c(-0.08, 0.08),
  ## range = c(-6.4, 6.4),
  ci_method = "HDI"
)
# Proportion of samples inside the ROPE [-0.08, 0.08]:

Parameter                   | inside ROPE
-----------------------------------------
scaleFERTILIZERscaleEQFALSE |      0.00 %
bayestestR::rope(fert.brm3,
  parameters = "FERTILIZER",
  range = c(-0.08, 0.08),
  ## range = c(-6.4, 6.4),
  ci_method = "HDI"
) |>
  plot()

bayestestR::equivalence_test(fert.brm3,
  parameters = "FERTILIZER",
  range = c(-0.08, 0.08),
  ## range = c(-6.4, 6.4),
  ci_method = "HDI"
)
# Test for Practical Equivalence

  ROPE: [-0.08 0.08]

Parameter                   |       H0 | inside ROPE |      95% HDI
-------------------------------------------------------------------
scaleFERTILIZERscaleEQFALSE | Rejected |      0.00 % | [0.60, 1.01]
bayestestR::equivalence_test(fert.brm3,
  parameters = "FERTILIZER",
  range = c(-0.08, 0.08),
  ## range = c(-6.4, 6.4),
  ci_method = "HDI"
) |>
  plot()
Picking joint bandwidth of 0.0159

newdata <- data.frame(FERTILIZER = c(200, 100))
fert.brm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  pairs()
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
 contrast                      estimate lower.HPD upper.HPD
 FERTILIZER200 - FERTILIZER100     80.6      59.9       100

Point estimate displayed: median 
HPD interval probability: 0.95 
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

Conclusions:

  • the change in Yield associated with an increase in Fertilizer concentration from 100 to 200 is 80.55.

We can alternatively express this as a percentage change…

newdata <- list(FERTILIZER = c(200, 100))
fert.brm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  regrid(transform = "log") |>
  pairs() |>
  regrid()
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
 contrast                    ratio lower.HPD upper.HPD
 FERTILIZER200/FERTILIZER100   1.6      1.41      1.82

Point estimate displayed: median 
HPD interval probability: 0.95 
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
# A tibble: 1 × 5
  variable                             median lower upper     P
  <chr>                                 <dbl> <dbl> <dbl> <dbl>
1 contrast FERTILIZER200/FERTILIZER100   1.60  1.41  1.82 0.880

If we want to derive other properties, such as the percentage change, then we use tidy_draws() and then simple tidyverse spreadsheet operation.

It is in combination with emmeans() that tidy_draws() really comes into its own.

fert.brm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  regrid(transform = "log") |>
  pairs() |>
  regrid() |>
  tidy_draws() |>
  summarise_draws(median,
    HDInterval::hdi,
    P = ~ mean(.x > 1.5)
  )
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
# A tibble: 1 × 5
  variable                             median lower upper     P
  <chr>                                 <dbl> <dbl> <dbl> <dbl>
1 contrast FERTILIZER200/FERTILIZER100   1.60  1.41  1.82 0.880
fert.mcmc <- fert.brm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  tidy_draws() |>
  rename_with(~ str_replace(., "FERTILIZER ", "p")) |>
  mutate(
    Eff = p200 - p100,
    PEff = 100 * Eff / p100
  )
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
fert.mcmc |> head()
# A tibble: 6 × 7
  .chain .iteration .draw  p200  p100   Eff  PEff
   <int>      <int> <int> <dbl> <dbl> <dbl> <dbl>
1      1          1     1  221.  132.  89.0  67.3
2      1          2     2  227.  131.  95.4  72.6
3      1          3     3  211.  131.  79.6  60.8
4      1          4     4  213.  131.  82.6  63.2
5      1          5     5  212.  134.  77.7  57.9
6      1          6     6  210.  132.  77.8  58.9

Now we can calculate the medians and HPD intervals of each column (and ignore the .chain, .iteration and .draw).

fert.mcmc |>
  tidyMCMC(
    estimate.method = "median",
    conf.int = TRUE, conf.method = "HPDinterval"
  )
# A tibble: 7 × 5
  term       estimate std.error conf.low conf.high
  <chr>         <dbl>     <dbl>    <dbl>     <dbl>
1 .chain          2       0.817      1         3  
2 .iteration    400.    231.         1       761  
3 .draw        1200.    693.         1      2281  
4 p200          214.      9.59     194.      232. 
5 p100          133.      8.38     117.      150. 
6 Eff            80.5    10.2       59.9     100. 
7 PEff           60.9    10.2       40.7      81.9

Alternatively, we could use median_hdci() to focus on a specific column.

fert.mcmc |> median_hdci(PEff)
# A tibble: 1 × 6
   PEff .lower .upper .width .point .interval
  <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1  60.1   40.7   81.9   0.95 median hdci     
fert.mcmc |> median_hdci(PEff, Eff)
# A tibble: 1 × 9
   PEff PEff.lower PEff.upper   Eff Eff.lower Eff.upper .width .point .interval
  <dbl>      <dbl>      <dbl> <dbl>     <dbl>     <dbl>  <dbl> <chr>  <chr>    
1  60.1       40.7       81.9  80.6      59.9      100.   0.95 median hdci     
fert.mcmc |>
  ggplot() +
  geom_density(aes(x = PEff))

Conclusions:

  • the percentage change in Yield associated with an increase in Fertilizer concentration from 100 to 200 is 60.15% with a 95% credibility interval of 40.66 81.94

To get the probability that the effect is greater than a 50% increase.

fert.mcmc |> summarise(P = mean(PEff > 50))
# A tibble: 1 × 1
      P
  <dbl>
1 0.880

Conclusions:

  • the probability that Yield will increase by more than 50 following an increase in Fertilizer concentration from 100 to 200 is 0.88.

Finally, we could alternatively use hypothesis(). Note that when we do so, the estimate is the difference between the effect and the hypothesised effect (50%).

fert.mcmc |> hypothesis("PEff>50")
Hypothesis Tests for class :
       Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
1 (PEff)-(50) > 0    10.94     10.24    -4.77    28.26        7.3      0.88
  Star
1     
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
newdata <- list(FERTILIZER = c(200, 100))
fert.brm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  pairs() |>
  tidy_draws() |>
  summarise(across(contains("contrast"),
    list(
      P = ~ mean(. > 80),
      HDCI = ~ median_hdci(.)
    ),
    .names = c("{.fn}")
  )) |>
  tidyr::unpack(HDCI)
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
# A tibble: 1 × 7
      P     y  ymin  ymax .width .point .interval
  <dbl> <dbl> <dbl> <dbl>  <dbl> <chr>  <chr>    
1 0.522  80.6  59.9  100.   0.95 median hdci     

We can also express this as a percentage change

newdata <- list(FERTILIZER = c(200, 100))
## Simple
fert.brm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  regrid(transform = "log") |>
  pairs() |>
  regrid()
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
 contrast                    ratio lower.HPD upper.HPD
 FERTILIZER200/FERTILIZER100   1.6      1.41      1.82

Point estimate displayed: median 
HPD interval probability: 0.95 
## More advanced (both P and percent change)
fert.mcmc <- fert.brm3 |>
  emmeans(~FERTILIZER, at = newdata) |>
  regrid(transform = "log") |>
  pairs() |>
  regrid() |>
  tidy_draws() |>
  mutate(across(contains("contrast"), ~ 100 * (. - 1)))
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
fert.mcmc |>
  summarise(across(contains("contrast"),
    list(
      P = ~ mean(. > 50),
      HDCI = ~ median_hdci(.)
    ),
    .names = c("{.fn}")
  )) |>
  tidyr::unpack(HDCI)
# A tibble: 1 × 7
      P     y  ymin  ymax .width .point .interval
  <dbl> <dbl> <dbl> <dbl>  <dbl> <chr>  <chr>    
1 0.880  60.1  40.7  81.9   0.95 median hdci     
fert.brm3 |>
  linpred_draws(newdata = as.data.frame(newdata)) |>
  ungroup() |>
  group_by(.draw) |>
  summarise(
    Eff = -1 * diff(.linpred),
    PEff = 100 * Eff / .linpred[2]
  ) |>
  ungroup() |>
  mutate(P = mean(PEff > 50)) |>
  pivot_longer(cols = -.draw) |>
  group_by(name) |>
  median_hdci()
# A tibble: 3 × 7
  name   value .lower  .upper .width .point .interval
  <chr>  <dbl>  <dbl>   <dbl>  <dbl> <chr>  <chr>    
1 Eff   80.6   59.9   100.      0.95 median hdci     
2 P      0.880  0.880   0.880   0.95 median hdci     
3 PEff  60.1   40.7    81.9     0.95 median hdci     
## OR
fert.brm3 |>
  epred_draws(newdata = as.data.frame(newdata)) |>
  ungroup() |>
  group_by(.draw) |>
  summarise(
    Eff = -1 * diff(.epred),
    PEff = 100 * Eff / .epred[2]
  ) |>
  ungroup() |>
  mutate(P = mean(PEff > 50)) |>
  pivot_longer(cols = -.draw) |>
  group_by(name) |>
  median_hdci()
# A tibble: 3 × 7
  name   value .lower  .upper .width .point .interval
  <chr>  <dbl>  <dbl>   <dbl>  <dbl> <chr>  <chr>    
1 Eff   80.6   59.9   100.      0.95 median hdci     
2 P      0.880  0.880   0.880   0.95 median hdci     
3 PEff  60.1   40.7    81.9     0.95 median hdci     
## OR for prediction of individual values
fert.brm3 |>
  predicted_draws(newdata = as.data.frame(newdata)) |>
  ungroup() |>
  group_by(.draw) |>
  summarise(
    Eff = -1 * diff(.prediction),
    PEff = 100 * Eff / .prediction[2]
  ) |>
  ungroup() |>
  mutate(P = mean(PEff > 50)) |>
  pivot_longer(cols = -.draw) |>
  group_by(name) |>
  median_hdci()
# A tibble: 3 × 7
  name   value .lower  .upper .width .point .interval
  <chr>  <dbl>  <dbl>   <dbl>  <dbl> <chr>  <chr>    
1 Eff   79.5    5.20  147.      0.95 median hdci     
2 P      0.627  0.627   0.627   0.95 median hdci     
3 PEff  59.6   -7.20  142.      0.95 median hdci     

12 Summary figures

fert.list <- with(fert, list(FERTILIZER = seq(min(FERTILIZER), max(FERTILIZER), len = 100)))
newdata <- emmeans(fert.rstanarm3, ~FERTILIZER, at = fert.list) |> as.data.frame()
head(newdata)
 FERTILIZER   emmean lower.HPD upper.HPD
   25.00000 72.18439  48.51879  102.5765
   27.27273 74.09538  50.54457  103.8918
   29.54545 75.97198  52.46619  104.9689
   31.81818 77.78670  54.64772  106.3487
   34.09091 79.62726  56.96628  107.8963
   36.36364 81.40451  59.28484  109.5472

Point estimate displayed: median 
HPD interval probability: 0.95 
ggplot(newdata, aes(y = emmean, x = FERTILIZER)) +
  geom_point(data = fert, aes(y = YIELD)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD), fill = "blue", alpha = 0.3) +
  scale_y_continuous("YIELD") +
  scale_x_continuous("FERTILIZER") +
  theme_classic()

## spaghetti plot
newdata <- emmeans(fert.rstanarm3, ~FERTILIZER, at = fert.list) |>
  gather_emmeans_draws()
newdata |> head()
# A tibble: 6 × 5
# Groups:   FERTILIZER [1]
  FERTILIZER .chain .iteration .draw .value
       <dbl>  <int>      <int> <int>  <dbl>
1         25     NA         NA     1   60.0
2         25     NA         NA     2   69.7
3         25     NA         NA     3   54.2
4         25     NA         NA     4   71.8
5         25     NA         NA     5   63.0
6         25     NA         NA     6   52.1
ggplot(newdata, aes(y = .value, x = FERTILIZER)) +
  geom_line(aes(group = .draw), alpha = 0.01) +
  geom_point(data = fert, aes(y = YIELD))

fert.list <- with(fert, list(FERTILIZER = seq(min(FERTILIZER), max(FERTILIZER), len = 100)))
newdata <- emmeans(fert.rstanarm3, ~FERTILIZER, at = fert.list) |> as.data.frame()
head(newdata)
 FERTILIZER   emmean lower.HPD upper.HPD
   25.00000 72.18439  48.51879  102.5765
   27.27273 74.09538  50.54457  103.8918
   29.54545 75.97198  52.46619  104.9689
   31.81818 77.78670  54.64772  106.3487
   34.09091 79.62726  56.96628  107.8963
   36.36364 81.40451  59.28484  109.5472

Point estimate displayed: median 
HPD interval probability: 0.95 
ggplot(newdata, aes(y = emmean, x = FERTILIZER)) +
  geom_point(data = fert, aes(y = YIELD)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD), fill = "blue", alpha = 0.3) +
  scale_y_continuous(expression(Grass ~ yield ~ (g.m^-2)), breaks = seq(50, 300, by = 50)) +
  scale_x_continuous(expression(Fertilizer ~ concentration ~ (g.m^-2))) +
  theme_classic()

## spaghetti plot
newdata <- emmeans(fert.rstanarm3, ~FERTILIZER, at = fert.list) |>
  gather_emmeans_draws()
newdata |> head()
# A tibble: 6 × 5
# Groups:   FERTILIZER [1]
  FERTILIZER .chain .iteration .draw .value
       <dbl>  <int>      <int> <int>  <dbl>
1         25     NA         NA     1   60.0
2         25     NA         NA     2   69.7
3         25     NA         NA     3   54.2
4         25     NA         NA     4   71.8
5         25     NA         NA     5   63.0
6         25     NA         NA     6   52.1
ggplot(newdata, aes(y = .value, x = FERTILIZER)) +
  geom_line(aes(group = .draw), colour = "orange", alpha = 0.01) +
  geom_point(data = fert, aes(y = YIELD)) +
  scale_y_continuous(expression(Grass ~ yield ~ (g.m^-2)), breaks = seq(50, 300, by = 50)) +
  scale_x_continuous(expression(Fertilizer ~ concentration ~ (g.m^-2))) +
  theme_classic()

fert.grid <- with(fert, list(FERTILIZER = modelr::seq_range(FERTILIZER, n = 100)))
fert.grid <- with(fert, list(FERTILIZER = seq(min(FERTILIZER), max(FERTILIZER), len = 100)))

newdata <- fert.brm3 |>
  emmeans(~FERTILIZER, at = fert.grid) |>
  as.data.frame()
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
head(newdata)
 FERTILIZER   emmean lower.HPD upper.HPD
   25.00000 73.11793  42.75258  98.64013
   27.27273 74.92301  45.71520 100.59468
   29.54545 76.76315  48.14779 102.16938
   31.81818 78.64609  50.36263 103.57037
   34.09091 80.46620  53.74772 106.39867
   36.36364 82.30375  54.49654 106.37236

Point estimate displayed: median 
HPD interval probability: 0.95 
ggplot(newdata, aes(y = emmean, x = FERTILIZER)) +
  geom_point(data = fert, aes(y = YIELD)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD), fill = "blue", alpha = 0.3) +
  scale_y_continuous("YIELD") +
  scale_x_continuous("FERTILIZER") +
  theme_classic()

## spaghetti plot
newdata <- emmeans(fert.brm3, ~FERTILIZER, at = fert.grid) |>
  gather_emmeans_draws()
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
newdata |> head()
# A tibble: 6 × 5
# Groups:   FERTILIZER [1]
  FERTILIZER .chain .iteration .draw .value
       <dbl>  <int>      <int> <int>  <dbl>
1         25     NA         NA     1   65.4
2         25     NA         NA     2   59.8
3         25     NA         NA     3   71.3
4         25     NA         NA     4   68.6
5         25     NA         NA     5   75.9
6         25     NA         NA     6   73.8
ggplot(newdata, aes(y = .value, x = FERTILIZER)) +
  geom_line(aes(group = .draw), color = "blue", alpha = 0.01) +
  geom_point(data = fert, aes(y = YIELD)) +
  theme_classic()

fert.list <- with(fert, list(FERTILIZER = seq(min(FERTILIZER), max(FERTILIZER), len = 100)))
newdata <- emmeans(fert.brm3, ~FERTILIZER, at = fert.list) |> as.data.frame()
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
head(newdata)
 FERTILIZER   emmean lower.HPD upper.HPD
   25.00000 73.11793  42.75258  98.64013
   27.27273 74.92301  45.71520 100.59468
   29.54545 76.76315  48.14779 102.16938
   31.81818 78.64609  50.36263 103.57037
   34.09091 80.46620  53.74772 106.39867
   36.36364 82.30375  54.49654 106.37236

Point estimate displayed: median 
HPD interval probability: 0.95 
ggplot(newdata, aes(y = emmean, x = FERTILIZER)) +
  geom_point(data = fert, aes(y = YIELD)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD), fill = "blue", alpha = 0.3) +
  scale_y_continuous(expression(Grass ~ yield ~ (g.m^-2)), breaks = seq(50, 300, by = 50)) +
  scale_x_continuous(expression(Fertilizer ~ concentration ~ (g.m^-2))) +
  theme_classic()

## spaghetti plot
newdata <- emmeans(fert.brm3, ~FERTILIZER, at = fert.list) |>
  gather_emmeans_draws()
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
newdata |> head()
# A tibble: 6 × 5
# Groups:   FERTILIZER [1]
  FERTILIZER .chain .iteration .draw .value
       <dbl>  <int>      <int> <int>  <dbl>
1         25     NA         NA     1   65.4
2         25     NA         NA     2   59.8
3         25     NA         NA     3   71.3
4         25     NA         NA     4   68.6
5         25     NA         NA     5   75.9
6         25     NA         NA     6   73.8
ggplot(newdata, aes(y = .value, x = FERTILIZER)) +
  geom_line(aes(group = .draw), colour = "orange", alpha = 0.01) +
  geom_point(data = fert, aes(y = YIELD)) +
  scale_y_continuous(expression(Grass ~ yield ~ (g.m^-2)), breaks = seq(50, 300, by = 50)) +
  scale_x_continuous(expression(Fertilizer ~ concentration ~ (g.m^-2))) +
  theme_classic()

13 Methods

  • The relationship between grass yield and fertiliser concentration was explored using linear regression in a Bayesian framework.
  • specifically, the yield of grass was modelled against fertiliser concentration with a Gaussian family and weakly informative priors.
  • (equation and priors)
  • the Bayesian model included 3 chains, each of 5000 iterations, thinned to a rate of 5, and excluded the first 1000 iterations (warmup)
  • the model was found to be well mixed and converged (all Rhat < 1.05) on a stable posterior (see suppl.) and was validated via DHARMa (cite) residuals
  • all statistical models were performed in the R (4.4.1) Statistical and Graphical Environment (R Core Team 2024) via the brms (Bürkner 2017) package

14 Results

  • Grass yield was found to have a positive linear relationship with fertiliser concentration (stats)
  • Every one unit increase in fertiliser concentration was found to be associated with a 0.8 gram increase in grass yield
  • Doubling the fertiliser concentration from 100 to 200 (units) is expected to increase the yield by x grams
fert.brm3 |> report()
Response residuals not available to calculate mean square error. (R)MSE
  is probably not reliable.
Start sampling
Response residuals not available to calculate mean square error. (R)MSE
  is probably not reliable.
We fitted a Bayesian linear model (estimated using MCMC sampling with 3 chains
of 5000 iterations and a warmup of 1000) to predict YIELD with FERTILIZER
(formula: YIELD ~ scale(FERTILIZER, scale = FALSE)). Priors over parameters
were set as normal (mean = 162.00, SD = 90.40) distributions. The model's
explanatory power is substantial (R2 = 0.92, 95% CI [0.82, 0.93], adj. R2 =
0.89).  Within this model:

  - The effect of b Intercept (Median = 163.41, 95% CI [148.73, 178.64]) has a
100.00% probability of being positive (> 0), 100.00% of being significant (>
3.20), and 100.00% of being large (> 19.19). The estimation successfully
converged (Rhat = 1.000) and the indices are reliable (ESS = 2456)
  - The effect of b scaleFERTILIZERscaleEQFALSE (Median = 0.81, 95% CI [0.60,
1.01]) has a 100.00% probability of being positive (> 0), 0.00% of being
significant (> 3.20), and 0.00% of being large (> 19.19). The estimation
successfully converged (Rhat = 1.000) and the indices are reliable (ESS = 2314)

Following the Sequential Effect eXistence and sIgnificance Testing (SEXIT)
framework, we report the median of the posterior distribution and its 95% CI
(Highest Density Interval), along the probability of direction (pd), the
probability of significance and the probability of being large. The thresholds
beyond which the effect is considered as significant (i.e., non-negligible) and
large are |3.20| and |19.19| (corresponding respectively to 0.05 and 0.30 of
the outcome's SD). Convergence and stability of the Bayesian sampling has been
assessed using R-hat, which should be below 1.01 (Vehtari et al., 2019), and
Effective Sample Size (ESS), which should be greater than 1000 (Burkner,
2017)., We fitted a Bayesian linear model (estimated using MCMC sampling with 3
chains of 5000 iterations and a warmup of 1000) to predict YIELD with
FERTILIZER (formula: YIELD ~ scale(FERTILIZER, scale = FALSE)). Priors over
parameters were set as normal (mean = 0.00, SD = 2.00) distributions. The
model's explanatory power is substantial (R2 = 0.92, 95% CI [0.82, 0.93], adj.
R2 = 0.89).  Within this model:

  - The effect of b Intercept (Median = 163.41, 95% CI [148.73, 178.64]) has a
100.00% probability of being positive (> 0), 100.00% of being significant (>
3.20), and 100.00% of being large (> 19.19). The estimation successfully
converged (Rhat = 1.000) and the indices are reliable (ESS = 2456)
  - The effect of b scaleFERTILIZERscaleEQFALSE (Median = 0.81, 95% CI [0.60,
1.01]) has a 100.00% probability of being positive (> 0), 0.00% of being
significant (> 3.20), and 0.00% of being large (> 19.19). The estimation
successfully converged (Rhat = 1.000) and the indices are reliable (ESS = 2314)

Following the Sequential Effect eXistence and sIgnificance Testing (SEXIT)
framework, we report the median of the posterior distribution and its 95% CI
(Highest Density Interval), along the probability of direction (pd), the
probability of significance and the probability of being large. The thresholds
beyond which the effect is considered as significant (i.e., non-negligible) and
large are |3.20| and |19.19| (corresponding respectively to 0.05 and 0.30 of
the outcome's SD). Convergence and stability of the Bayesian sampling has been
assessed using R-hat, which should be below 1.01 (Vehtari et al., 2019), and
Effective Sample Size (ESS), which should be greater than 1000 (Burkner, 2017).
and We fitted a Bayesian linear model (estimated using MCMC sampling with 3
chains of 5000 iterations and a warmup of 1000) to predict YIELD with
FERTILIZER (formula: YIELD ~ scale(FERTILIZER, scale = FALSE)). Priors over
parameters were set as student_t (location = 0.00, scale = 90.40)
distributions. The model's explanatory power is substantial (R2 = 0.92, 95% CI
[0.82, 0.93], adj. R2 = 0.89).  Within this model:

  - The effect of b Intercept (Median = 163.41, 95% CI [148.73, 178.64]) has a
100.00% probability of being positive (> 0), 100.00% of being significant (>
3.20), and 100.00% of being large (> 19.19). The estimation successfully
converged (Rhat = 1.000) and the indices are reliable (ESS = 2456)
  - The effect of b scaleFERTILIZERscaleEQFALSE (Median = 0.81, 95% CI [0.60,
1.01]) has a 100.00% probability of being positive (> 0), 0.00% of being
significant (> 3.20), and 0.00% of being large (> 19.19). The estimation
successfully converged (Rhat = 1.000) and the indices are reliable (ESS = 2314)

Following the Sequential Effect eXistence and sIgnificance Testing (SEXIT)
framework, we report the median of the posterior distribution and its 95% CI
(Highest Density Interval), along the probability of direction (pd), the
probability of significance and the probability of being large. The thresholds
beyond which the effect is considered as significant (i.e., non-negligible) and
large are |3.20| and |19.19| (corresponding respectively to 0.05 and 0.30 of
the outcome's SD). Convergence and stability of the Bayesian sampling has been
assessed using R-hat, which should be below 1.01 (Vehtari et al., 2019), and
Effective Sample Size (ESS), which should be greater than 1000 (Burkner, 2017).

15 References

Bürkner, Paul-Christian. 2017. brms: An R Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software 80 (1): 1–28. https://doi.org/10.18637/jss.v080.i01.
Cohen, J. 1988. Statistical Power Analysis for the Behavioral Sciences. Lawrence Erlbaum Associates.
Fowler, J., L. Cohen, and P. Jarvis. 1998. Practical Statistics for Field Biology. England: John Wiley & Sons.
Kruschke, John K. 2018. “Rejecting or Accepting Parameter Values in Bayesian Estimation.” Advances in Methods and Practices in Psychological Science 1 (2): 270–80. https://doi.org/10.1177/2515245918771304.
R Core Team. 2024. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.